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1.  PROJECT  SUMMARY 


The  overall  goal  of  the  work  conducted  in  this  project  was  to  develop  a  combined  simulation  and 
experimental  capability  to  further  the  understanding  of  unsteady  wing  aerodynamics  of  flapping  micro  air 
vehicles  (MAVs)  and  to  assist  in  design  of  such  vehicles*  The  primary  objective  of  the  overall  effort  was  to 
develop  and  validate  a  simulation  environment  with  high-fidelity  modeling  of  flapping-wing  MAV  flight 
physics.  The  next  objective  was  to  study  the  fluid  physics  associated  with  low  Reynolds  number  flapping 
wings  to  facilitate  the  design  of  effective  high-bandwidth  control  of  wing  beat  kinematics  for  MAVs* 
Another  objective  was  to  build  a  fundamental  knowledge  base  for  the  physics  of  flapping  wing  Micro  Air 
Vehicles  (MAVs)*  The  Phase  2  work  had  two  major  components: 

•  Simulation  environment  based  on  an  advanced  Computational  Fluid  Dynamics  (CFD)  solver  called 
Loci-STREAM,  and 

*  Experimental  effort  to  validate  the  simulations  and  to  guide  model  development. 


Loci -STREAM  served  as  the  simulation  tool  for  investigating  flapping  wing  aerodynamics.  Several  pitching 
and  plunging  airfoil  cases  were  investigated  during  this  project  which  demonstrated  the  need  for  accurate  and 
efficient  CFD  capability  to  resolve  the  unsteady  How  fields  associated  with  flapping  wings.  The  simulations 
demonstrated  that  a  variation  of  the  Reynolds  number  (wing  sizing,  flapping  frequency,  etc)  leads  to  changes 
in  the  leading  edge  vortex  (LEV)  and  span  wise  flow  structures,  which  impact  the  aerodynamic  force 
generation*  While  in  classical  stationary  wing  theory  the  tip  vortices  (TiVs)  are  seen  as  wasted  energy,  in 
flapping  flight,  they  can  interact  with  the  LEV  to  enhance  lift  without  increasing  the  power  requirements* 
Surrogate  modeling  techniques  can  assess  the  aerodynamic  outcomes  between  two-  and  three-dimensional 
wings.  The  combined  effect  of  the  TiVs,  the  LEV,  and  jet  can  improve  the  aerodynamics  of  a  flapping  wing* 
Efforts  are  continuing  in  broadening  Loci-STREAM 's  capability  to  address  aeroelasticity.  In  particular, 
chord  wise  flexibility  in  the  forward  flight  can  substantially  adjust  the  projected  area  normal  to  the  flight 
trajectory  via  shape  deformation,  hence  redistributing  thrust  and  lift*  Span  wise  flexibility  in  the  forward 
flight  creates  shape  deformation  from  the  wing  root  to  the  wing  tip  resulting  in  varied  phase  shift  and 
effective  angle  of  attack  distribution  along  the  wing  span*  Numerous  open  issues  in  flapping  wing 
aerodynamics  along  these  lines  exist.  Additionally,  numerous  issues  exist  with  regards  to  the  computational 
techniques  needed  to  model  such  flow  fields.  Several  such  issues  were  addressed  during  this  project  such  as 
advanced  grid  movement  techniques,  fluid-structure  interaction  capability,  etc*  While  the  simulation  tool  is 
still  being  developed  and  refined,  significant  advancement  is  expected  once  the  tool  reaches  the  intended 
potential  in  the  near  future. 


On  the  experimental  side,  investigations  were  conducted  on  model  geometries  to  validate  the  CFD  capability 
and  to  assist  in  model  development.  Pitching  and  plunging  airfoil  at  Reynolds  number,  *Re=  10,000,  and  pure 
sinusoidal  effective  angle  of  attack  motion  has  been  investigated.  The  experiments  were  conducted  at  the 
University  of  Michigan  low-turbulence  water  channel  facility  using  2D  phase-averaged  particle  image 
velocimetry  (PIV).  The  effect  of  non-dimensional  parameters  governing  pitching  and  plunging  motion  such 
as  Strouhal  number  (St),  reduced  frequency  (ft),  and  the  plunge  amplitude  (ho)  was  investigated  for  the  same 
effective  angle  of  attack  kinematics*  The  formation  of  a  leading  edge  vortex  (LEV)  and  a  trailing  edge  vortex 
(TEV)  was  observed  for  all  the  cases  studied*  The  formation  phase  of  the  LEV  was  found  to  be  dependent  on 
k ;  the  LEV  formation  is  delayed  for  higher  k  value*  It  was  found  that  for  cases  with  the  same  k  the  velocity 
profiles  normal  to  the  airfoil  surface  closely  follow  each  other  in  all  cases  independent  of  pitch  rate  and  pivot 
point  effect.  Analysis  on  the  locations  of  the  LEV  core  based  on  the  Q-criterion  and  local  streamline  patterns 
helped  identify  the  trajectory  of  the  LEV  core  with  respect  to  the  airfoil.  Additionally,  a  trend  in  the  LEV 
circulation  was  observed* 
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2,  DEVELOPMENT  OF  COMPUTATIONAL  FRAMEWORK 


2.1  A  High-Fidelity  Parallel  Simulation  Environment  for  Flapping  Wing  MA  Vs: 
Loci-STREAM 

One  of  the  main  objectives  of  this  project  was  to  demonstrate  an  efficient,  high-fidelity  simulation  tool  which 
can  be  used  to  enhance  the  understanding  of  MAV  fluid  physics  and  to  assist  in  the  design  of  such  vehicles. 
A  high-fidelity  simulation  capability  requires  the  use  of  a  full  unsteady,  viscous,  Navier-Stokes  based  CFD 
methodology.  Over  the  past  decade,  such  a  methodology  has  been  developed  by  Streamline  Numerics 
resulting  in  a  code  called  Loci-STREAM  for  large  scale  parallel  CFD  applications  (with  funding  from 
NASA).  The  Loci-STREAM  code  which  has  been  developed  within  a  novel  programming  framework  called 
Loci  (developed  at  the  Mississippi  State  University).  Loci-STREAM  has  several  key  features  which  make  it 
a  very  attractive  framework  as  an  MAV  flow  simulation  and  design  tool: 

*  The  basic  algorithm  is  an  unsteady,  pressure-based  finite  volume  method,  which  is  ideal  for  low 
Mach  number  flows  including  truly  incompressible  flows,  which  are  typical  of  M  A  Vs. 

*  The  algorithm  has  been  developed  for  generalized  unstructured  grids  (employing  arbitrarily  shaped 
polyhedra)  which  allows  great  flexibility  and  economy  in  grid  generation  for  complex  geometries 
such  as  MAVs  (involving  wing-body  junctions). 

*  The  Loci  programming  environment  facilitates  efficient  large-scale  parallel  computing  which  will 
allow  the  simulation  of  grid-independent,  three-dimensional  unsteady  flows  around  MAVs  with  fast 
turnaround  times. 

The  work  involved  developing  the  foundation  for  the  advanced  high-fidelity  simulation  environment  in  the 
Loci-STREAM  CFD  software  package,  involving  advanced  mesh  movement  and  robust  unsteady  modeling 
capabilities  that  are  essential  for  modeling  MAV  flight  physics.  The  basic  grid  movement  algorithm  was 
incorporated  into  Loci-STREAM.  Following  that,  test  cases  for  the  translation  and  rotation  of  both  thin 
membranes  and  airfoils  were  used  to  demonstrate  that  the  grid  movement  algorithm  is  capable  of  allowing 
large  geometric  deflections  (on  the  order  of  the  body  size)  even  in  the  presence  of  highly-refined  viscous 
grids,  while  still  maintaining  sufficient  grid  quality. 

2.2  Loci-STREAM:  Algorithm  for  Moving  Boundary  Problems 

Loci-STREAM  is  a  parallelized  unstructured  curvilinear  pressure- based  finite- volume  code  with  moving  grid 
capabilities.  In  the  algorithm  employed  in  Loci-STREAM,  the  momentum  and  pressure  equations  are 
segregated  and  the  solution  is  advanced  using  an  implicit  geometrically  conservative  time  integration 
method.  The  convection  terms  are  treated  using  the  second  order  upwind  scheme  while  pressure  and  viscous 
terms  are  treated  using  second  order  schemes.  The  geometric  conservation  law,  a  necessary  consideration  in 
domains  with  moving  boundaries,  is  satisfied.A  detailed  description  of  the  method  implemented  in  two- 
dimensions  is  given  elsewhere  (Wright  and  Smith  2001,  Smith  and  Wright  2004).  Here  we  present  the  main 
governing  equations  very  briefly;  the  governing  equations  for  the  numerical  simulation  are  the  RANS 
equations  coupled  with  Menter’s  SST  model  (Menter  1994),  and  the  continuity  equation  for  incompressible 
flow. 
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where  altp9fJ*9  y,  rrk,  o^,F2  are  defined  as  in  M enter’s  SST  formulation  is  the  velocity  component  in  the 
itU  direction,  Xj  is  the  ith  component  of  the  position  vector,  i  is  time,  p  is  density,  p  is  pressure,  v  is  the 
kinematic  viscosity,  vt  is  the  eddy  viscosity,  S  —  ^fZS^S^j  is  the  invariant  measure  of  the  strain  rate* 
Compared  to  Menter’s  original  SST  turbulence  model  a  limiter  has  been  built  in  (Menter  2003)  to  the 
production  term,  Pk,  in  the  turbulence  kinetic  energy  (TKE)  equation,  Eq.  (1),  as 


du:  fdu ;  duA 

Pk  =  fit’dx^\d7}  +  dxl)‘ 

Pk  —  min(Pkl  10  ■  {3*pko))t 


where  Pk  is  the  production  term  in  the  original  SST  formulation,  to  prevent  the  build-up  of  turbulence  in 
stagnation  regions.  Another  change  is  the  use  of  invariant  measure  of  the  strain-rate  tensor  in  the  formulation 
for  the  eddy  viscosity  instead  of  the  vorticity  magnitude,  H  =  The  strain-rate  invariant  is 

considered  to  be  a  better  measure  for  the  fluid  deformation,  since  the  Boussinesq  approximation  is  also  based 
on  the  strain-rate.  The  two  differences  between  the  original  and  the  modified  SST  formulation  are 
summarized  in  Table  2.2-1 . 
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Vt  — - 

max(a1w,flF2) 

a^k 

y .  — - 

maxCd,  <u,S  F2) 

Loci-STREAM  has  been  demonstrated  to  exhibit  near-linear  scalability  for  general 
incompressible/compressible,  turbulent  flows  using  unstructured  meshes.  Loci-STREAM  operates  with 
near- linear  scalability  for  both  shared-memory  and  distributed-memory  computer  architectures,  as  shown  in 
Figure  2.2-1. 
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Figure  2.2-1.  Speed-up  ratio  for  a  typical  problem  using  Loci-STREAM  on:  (a)  shared  memory  architecture, 
and  (b)  distributed  memory  architecture. 


The  equation  governing  the  static  equilibrium  of  a  homogeneous  linear  elastic  solid  in  the  absence  of  body 
forces  expressed  in  terms  of  displacements  is  the  Navier-Cauchy  equation 

(1-— )V(V-w)  +  V2m  =  0  (2.2) 

2  v 

where  u  is  the  small  displacement  vector  measured  with  respect  to  an  undeformed  reference  configuration 
and  vis  Poisson’s  ratio.  Solutions  to  Eq.  (2.2),  obtained  under  appropriate  boundary  conditions,  are  used  to 
update  the  nodal  coordinates  according  to 

at"*1  =  x”  +  u  (2.3) 

where  1  and  x"  are  the  nodal  coordinates  at  the  current  and  previous  time  levels,  respectively.  For  portions 
of  the  computational  domain  boundary  subject  to  motion,  Dirichlet  boundary  conditions  are  imposed 
according  to 

«  =  g(0  on  tf  (2-4) 

where  g(t)  prescribes  the  displacement  of  boundary  nodes  and  P/  is  the  moving  portion  of  the  domain 
boundary.  Other  boundary  conditions  such  as  stress-free  or  symmetry  conditions  may  be  imposed  on  other 
portions  of  the  domain  boundary  to  accommodate  specific  kinematic  and  geometric  situations. 


Integrating  Eq.  (2.2)  over  a  the  same  control  volume  used  to  discretize  the  momentum  equation  gives  the 
Navier-Cauchy  equation  in  weak  form  as 


1 

(l-2u) 


«jV(V  •  u)dQ  +  c[(V  ■  Vh  )dCl  =  0 

ad  a,/ 


(2.5) 


where  is  the  control  volume. 


Following  Stein  and  Tezduyar  (2002),  the  stiffness  of  an  individual  control  volume  in  the  mesh  is  modified 
based  on  its  volume  at  the  current  time  step.  In  a  face/edge-based  finite  volume  method  such  as  the  one  used 
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here,  control  volumes  are  modified  on  an  edge  basis  according  to 


where  x  is  the  stiffening  parameter,  is  the  edge  contribution  to  the  dual  control  volume  and  Hrvi  is  an 
arbitrary  reference  volume  introduced  to  maintain  dimensional  consistency.  When  %  is  zero,  Eqs.  (6)  and 
(13)  reduce  to  the  system  governing  the  static  response  of  a  homogeneous  linear  elastic  solid  to  imposed 
displacement  boundary  conditions.  When  %  takes  on  values  greater  than  zero  the  discrete  system  mimics  the 
behavior  of  a  non-homogenous  linear  elastic  system  where  the  local  material  stiffness  has  been  increased  by 
a  factor  (7/fly*  The  overall  effect  of  this  volume- based  stiffening  is  to  cause  small  control  volumes  near  a 
moving  boundary  to  exhibit  approximately  rigid  body  motion  while  exporting  the  bulk  of  the  deformation  to 
regions  with  larger  cell  sizes  where  deformation  can  be  more  readily  tolerated  without  severely  degrading  the 
mesh  quality. 

Coding  was  added  to  Loci-STREAM  to  satisfy  the  so-called  Geometric  Conservation  Law  (GCL)  which  is 
required  when  solving  fluid  dynamics  problems  on  generalized  moving  meshes.  In  addition,  a  general 
capability  has  also  been  implemented  to  specify  rigid-body  translation  and  rotation  grid  motions  required  for 
problems  of  interest  in  this  STTR  (plunging  and  pitching  airfoils)  project,  in  order  to  test  the  new  capability, 
the  flow  about  a  cylinder  of  diameter  D,  oscillating  in  a  stationary  fluid  has  been  used  as  a  benchmark 
(Uzunoglu  et  aL  2001).  The  transverse  motion  of  the  cylinder  is  given  by 

*(0  =  A  sin(Q/),  *(/)  =  /tQcos(Q0  “  UM  cos(Q/)  (2,7) 


Figure  2,2-2,  litviscid  force  on  oscillating  cyinder. 


where  A  is  the  amplitude  of  the  motion  and  fl  is  the  circular  frequency.  For  inviscid  flow  of  small  amplitude, 
the  force  per  unit  length  acting  on  the  cylinder  is  given  as  follows: 

F{t)=  Q.25p7rD*x(t)  (2.8) 

A  numerical  solution  was  obtained  with  Loci-STREAM  using  a  61x27  grid  An  amplitude,  A=0.l><D,  was 
used  which  is  considered  sufficiently  small  to  satisfy  the  conditions  of  the  analytical  benchmark.  Figure 
2.2-2  shows  a  comparison  of  the  analytical  and  numerical  solutions  for  ten  cycles.  It  is  clear  that  the 
numerical  solution  matches  the  benchmark  extremely  well. 
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23  Incorporation  of  Advanced  Grid  Movement  Techniques 

Phase  I  work  demonstrated  the  feasibility  of  Loci-STREAM  in  serving  as  a  simulation  tool  for  flapping  wing 
aerodynamics.  Several  2-D  pitching  and  plunging  airfoil  cases  were  investigated  during  Phase  1  work  to 
demonstrate  the  need  for  accurate  and  efficient  CFD  capability  to  resolve  the  unsteady  flow  fields  associated 
with  flapping  wings.  To  robustly  and  efficiently  simulate  the  aerodynamics  of  flapping-wing  MA Vs,  a  robust 
grid  movement  strategy  -  which  will  not  fait  under  conditions  of  extreme  grid  movement  -  is  a  key 
ingredient. 

2.3-1  Rigid-Body  Source  Term  Strategy  to  Reduce  Grid  Skewness  During  Movement 


In  our  previous  work  (during  Phase  t)  it  was  observed  that  the  original  Stein  and  Tezduyar  (2002)  grid 
movement  algorithm  (Smith  and  Wright  2005),  failed  to  produce  the  desired  results  when  the  so-called 
Jacobian-based  stiffening  was  activated.  Specifically,  tor  cases  involving  rigid-body  rotation  of  a  body 


Figure  2-3-1.  Results  for  three  test  cases  (Stein  and  Tezduyar  2002)  with  the  modified  grid  movement 
algorithm- 

surface,  grid  lines  in  the  vicinity  of  the  surface  failed  to  remain  orthogonal  to  the  surface,  resulting  in  skewed 
grid  elements  near  the  surface.  For  large  angles  of  rotation,  correspondingly  large  skewness  was  observed, 
rendering  the  mesh  essentially  useless  for  accurate  CFD  calculations. 

One  strategy  to  counteract  this  deficiency  is  to  use  explicit  rigid-body  rotation  for  the  grid  nodes  near  a 
rotating  surface.  While  this  strategy  is  acceptable  for  simple  two-dimensional  problems,  logistical  problems 
arise  for  general  three-dimensional  cases  where  different  parts  of  the  surface  (a  flapping  bird  for  example) 
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are  moving  in  different  ways,  and  one  must  identify  near-surface  nodes  with  a  unique  rigid-body  motion.  At 
places  in  the  mesh  such  as  wing-body  intersections,  this  strategy  appears  to  become  completely  intractable. 

In  lieu  of  the  above  approach,  a  simple  heuristic  strategy  of  adding  a  rigid-body  source  term  to  the  original 
linear  elasticity  equations  governing  mesh  movement  was  investigated.  The  source  term  takes  the  form  of  a 
diffusion  flux  of  a  modified  displacement  gradient  tensor,  in  which  the  diagonal  entries  have  been  zeroed 
out,  and  the  off-diagonal  entries  include  the  strain-free  compatibility  conditions.  This  modified  source  term 
approach  is  completely  general,  allowing  its  applicability  to  both  two-dimensional  and  three-  dimensional 
configurations  with  no  algorithmic  complexity.  Figure  23-1  shows  the  results  for  the  three  basic  Stein  and 
Tezduyar  (2002)  test  cases  obtained  with  this  modified  algorithm.  The  results  are  essentially  identical  to  the 
benchmark  calculations  and  show  a  great  improvement  in  the  ability  of  the  algorithm  to  maintain 
orthogonality  of  the  grid  lines  in  the  vicinity  of  a  rotating  surface. 


To  test  the  new'  algorithm  on  a  realistic  two-dimensional  configuration,  a  NACA0012  airfoil  grid  was 
generated,  suitable  for  a  high  Reynolds  number  turbulence  flow  problem,  with  wall  spacing  of  0.0001. 
Figure  23-2  shows  the  deformed  mesh  after  the  airfoil  has  moved  through  45  degrees  of  rotation.  Only 
boundary  node  displacements  have  been  assigned,  while  all  internal  node  displacements  are  computed  using 
the  new  mesh  movement  algorithm.  From  the  far  view,  one  cannot  tell  that  the  mesh  was  not  explicitly 
rotated.  The  near  view  shows  that  mesh  quality  and  near-orthogonality  have  been  maintained  even  for  this 
fine  spacing  near  the  surface  of  the  airfoil.  Overall,  the  new  modified  mesh  movement  algorithm  appears  to 
produce  high-quality  meshes  even  under  severe  boundary  deformations,  with  robustness  independent  of  the 
underlying  geometry. 


Figure  23-2.  Deformed  mesh  for  the  NACA  0012  airfoil  using  the  new  mesh  movement  algorithm. 


The  new  mesh  movement  strategy  in  Loci-STREAM  was  employed  for  flapping  wing  geometries  presented 
in  the  following  sections  of  this  report. 
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2,3,2  Strategies  for  Grid  Movement  Boundary  Condition  Specification 

One  of  the  primary  goals  of  the  present  effort  was  to  validate  Loci-STREAM  with  experimental  data  for 
flapping  wing  configurations.  The  experiments  conducted  at  AFRL  with  the  SD7003  airfoil  employed  a  3-D 
wing  inside  a  water  tunnel,  with  small  gaps  (on  the  order  of  1  mm)  between  the  wing  tips  and  the  water  tunnel 
walls  in  order  to  reduce  3-D  effects.  Since  3-D  effects  cannot  be  entirely  eliminated  in  the  water  tunnel,  it  is 
important  to  quantify  these  effects  so  that  one  may  understand  whether  deviations  from  experiment  are 
caused  by  numericai/computational  modeling  issues  or  by  differences  in  the  experimental  measurements  that 
are  actually  due  to  3-D  effects.  Thus,  detailed  simulations  of  the  complete  3-D  SD7003  wing  in  the  water 
tunnel  are  required.  Such  computations  are  very  challenging  due  to  the  presence  of  the  minute  gaps  between 
the  wing-tips  and  the  tunnel  walls,  and  a  robust  grid-movement  strategy  is  essential  in  order  to  maintain  grid 


Figure  2.3-3.  Global  view  of  the  computational  mesh  for  the  2-D  test  problem  which  consists  of  a  square  of 
the  same  dimension  as  the  water  tunnel  cross-section  (0,61m  x  0.61m)  and  a  finite  thickness  plate  (0,60in  x 
0.01m)  which  undergoes  a  periodic  oscillation. 

quality  under  the  presence  of  relatively  large  wing  displacements  compared  to  gap  dimensions.  First,  a  2-D 
model  problem  simulation  was  undertaken  in  order  to: 
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(a)  Determine  if  the  grid  generation  software  SolidMesh  (developed  at  Mississippi  State  University) 
currently  being  used  is  capable  of  generating  a  high-quality  mesh  in  the  gap  region  which  maintains 
appropriate  boundary  layer  resolution  while  providing  the  proper  isotropic  mesh  element  fill  between 
the  tunnel  wall  and  wing-tip  boundary  layers. 

(b)  Determine  the  best  strategies  for  grid  movement  given  the  presence  of  small  length-scale  gaps,  and 

(c)  Investigate  if  the  basic  grid-movement  algorithm  which  has  been  developed,  in  conjunction  with 
these  strategies  will  allow  us  to  maintain  grid  quality  in  the  gap  region  over  a  period  of  many  pitch 
cycles. 

Figure  2,3-3  hows  a  global  view  of  the  computational  mesh  for  the  2-D  test  problem  which  consists  of  a 
square  of  the  same  dimension  as  the  water  tunnel  cross-section  (0.61m  x  0.61m)  and  a  finite  thickness  plate 
(0.60m  x  0,01m)  which  undergoes  a  periodic  oscillation.  The  gap  between  the  plate  and  the  walls  of  The 
square  is  0.001m  which  corresponds  exactly  to  the  real  3-D  simulation.  Figure  2.3-4  shows  a  blow-up  of  the 
gap  region  where  one  can  see  boundary  layer  meshes  with  initial  spacing  of  1.0e~05m  growing  off  both  the 
plate  (wing  tip)  and  square  (tunnel  wall)  with  isotropic  cell  fill-in  where  the  boundary  layers  end.  Other  than 
a  small  region  of  the  wall  boundary  layer  where  quadrilateral  cells  have  been  collapsed  into  triangles  (no 
combination  of  grid-generation  parameters  would  eliminate  this),  SolidMesh  is  clearly  capable  of  generating 
the  required  mesh  in  this  region. 

In  the  experimental  cases  of  McGowan  et  ak  (2008)  a  wing  of  chord  0,1524m  was  put  through  a  series  of 
pure  plunges  with  an  amplitude  of  5%  of  the  chord  length  (0.00762m).  The  approximate  maximum  thickness 
of  the  SD7003  airfoil  based  on  this  chord  length  is  0.012m,  which  motivated  the  choice  of  a  plate  thickness 
of  0.01  m.  A  numerical  simulation  was  performed  in  which  the  plate  undergoes  an  oscillation  about  y=0.0  of 
magnitude  0.0 1  m  (approximately  30%  larger  than  actual  experiment).  Two  strategies  for  specifying  the 
boundary  conditions  of  the  grid  movement  were  examined,  and  are  detailed  as  follows: 


Figure  23-4,  Blow-up  of  the  gap  region  showing  boundary  layer  meshes  with  initial  spacing  of  1.0e-05m 
growing  off  both  the  plate  (wing  tip)  and  square  (tunnel  wall)  with  isotropic  cell  fill-in  where  the  boundary 
layers  end. 
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1)  The  motion  of  grid  nodes  on  the  plate  is  specified  with  Dirichlet  conditions  and  all  other  nodes  are 
computed  using  the  mesh  movement  algorithm.  All  wall  nodes  are  constrained  to  stay  in  the  plane 
of  the  wall. 

2)  Same  condition  for  plate  nodes.  Boundary  nodes  on  a  small  section  of  the  tunnel  wall  are  also 
prescribed  to  move  in  unison  with  the  airfoil.  Ail  other  wall  nodes  are  computed  using  the  mesh 
movement  algorithm  and  constrained  to  stay  in  the  plane  of  the  wall. 


Upon  investigation  of  strategy  I,  it  quickly  became  apparent  that  the  grid-movement  algorithm  was  not 
capable  of  dragging  the  nodes  along  the  wall  adjacent  to  the  airfoil  in  a  manner  sufficient  to  maintain 
adequate  grid  quality.  This  strategy  was  therefore  abandoned.  Using  strategy  2,  in  which  the  adjacent  wall 
nodes  are  specified  to  mirror  the  motion  of  the  airfoil,  it  was  found  that  all  nodes  between  the  wing  tip  and 
the  tunnel  wall  move  with  an  effective  rigid-body  motion  which  maintains  the  initial  quality  of  the  mesh.  The 
simulation  using  strategy  2  was  carried  out  for  a  total  duration  of  four  cycles  at  which  it  was  terminated  since 
it  became  readily  apparent  that  no  significant  accumulative  degradation  in  grid  quality  was  occurring.  Figure 
2.3-5  shows  a  close-up  of  the  upper  part  of  the  gap  region  when  the  plate  is  at  its  maximum  positive 
displacement  (y=O.0 1  m)  during  the  fourth  cycle.  Figure  2.3-6  shows  a  similar  shot  of  the  lower  part  the  gap 
region  when  the  plate  is  at  its  maximum  negative  displacement  (y=-0.01m)  during  the  fourth  cycle.  It  is  clear 
that  the  current  grid  movement  algorithm  (Stein  and  Tezduyar  2002)  along  with  strategy  2  is  very  effective  in 
maintaining  excellent  grid  quality  under  amplitudes  of  motion  of  current  interest.  We  also  anticipate  that 
these  results  will  directly  extend  to  the  full  3-D  wing  simulation  at  comparable  and  even  higher  amplitudes  of 
wing  motion. 


Figure  2.3-5.  Close-up  of  the  upper  part  of  the  gap  region  when  the  plate  is  at  its  maximum  positive 
displacement  (y=0.01  m)  during  the  fourth  cycle. 
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X 

Figure  2*3-6,  Lower  part  the  gap  region  when  the  plate  is  at  its  maximum  negative  displacement  (y=-G.0lm) 
during  the  fourth  cycle. 


23,3  Investigation  of  Stiffness  Parameter 

An  important  part  of  this  project  in  the  understanding  of  flapping- wing  flight  is  the  use  of  both  CFD  and 
experiment  in  a  symbiotic  fashion  to  investigate  the  fluid-dynamic  mechanisms  of  interest.  Initially,  the 
CFD  simulations  in  this  project  were  performed  for  unbounded  domains,  with  comparisons  made  to 
experimental  results  which  have  been  conducted  in  bounded  domains.  The  next  step  in  the  process  was  to 
assess  the  impact  of  bounded  domains  on  the  experimental  results  and  determine  if  there  are  any  significant 
effects  which  can  lead  to  deviations  between  CFD  and  experiment. 

As  a  starting  point  in  this  process,  a  grid  was  generated  for  RTO  case  (lb)  in  the  bounded  domain  of  a  water 
tunnel  with  no-slip  walls  both  above  and  below  the  airfoil.  While  the  top  boundary  is  in  reality  a  free  surface, 
the  use  of  no-slip  boundaries  above  and  below  simulates  the  case  in  which  plates  have  been  placed  on  the  top 
of  the  water  tunnel  to  suppress  the  effects  of  the  free  surface.  At  a  later  time,  the  actual  tree-surface  problem 
will  be  included  in  the  CFD  simulation.  Figure  2.3-7  shows  a  portion  of  the  computational  domain  around 
the  airfoil  near  its  initial  configuration.  The  overall  grid  extends  from  an  upstream  position  of  x=-2.0  to  a 
downstream  position  of  x=3.0,  while  the  lower  wall  lies  at  position  y=- 0,305  and  the  upper  wall  at  position 
y=0.305.  The  airfoil  chord  is  0.1524,  and  the  sinusoidal  motion  takes  the  airfoil  a  half  chord  length  about  the 
mean  position  >  =0,0. 
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Initial  investigations  have  concentrated  on  testing  the  grid  movement  algorithm  over  the  complete  range  of 
motion  of  the  airfoil  to  determine  the  grid  movement  parameters  required  for  maintaining  a  quality  mesh 
over  the  duration  of  the  simulation*  The  primary  parameter  of  interest  in  the  Stein  and  Tezduyar  (2002)  grid 


Figure  2.3-7. .  Computational  domain  around  the  airfoil  near  its  initial  configuration. 


movement  algorithm  is  the  stiffness  parameter,  which  scales  the  size  of  the  differential  stiffness  between 
small  and  large  grid  elements.  The  fundamental  idea  behind  the  algorithm  is  to  enhance  the  stiffness  of 
smaller  grid  elements  (which  tend  to  be  near  the  moving  body)  and  to  relax  the  stiffness  of  larger  elements 
(which  tend  to  be  away  from  the  moving  body). 

For  moving  mesh  cases  with  unbounded  domains,  relatively  large  values  of  the  stiffness  parameter  x  can  be 
used.  This  results  in  a  mesh  which  is  nearly  rigid  in  the  vicinity  of  the  moving  body,  while  allowing 
significant  mesh  deformation  to  occur  far  away  from  the  body  in  order  to  accommodate  large  physical 
displacements  of  the  moving  body  while  suffering  virtually  no  loss  in  grid  quality  at  any  location  in  the 
domain.  For  the  current  bounded  domain  case,  initial  investigations  have  shown  that  much  smaller  values  of 
the  stiffness  parameter  must  be  used  in  order  to  maintain  adequate  grid  quality. 

Figure  2*3-8  (a,b,c)  show  portions  of  the  CFD  mesh  for  x=F5  when  the  airfoil  is  located  at  its  maximum 
vertical  displacement*  For  this  relatively  large  value  of  the  stiffness  parameter,  the  grid  quality  in  the  vicinity 
of  the  leading  and  trailing  edges  of  the  airfoil  is  seen  to  be  excellent,  virtually  the  same  quality  as  the  initial 
configuration.  Mear  the  lower  wall  however,  the  mesh  has  undergone  what  appears  to  be  a  delamination 
process,  resulting  in  an  unacceptable  mesh  at  this  point  in  the  airfoil  displacement  cycle.  This  breakdown  in 
the  mesh  is  essentially  caused  by  too  high  a  stiffness  parameter,  in  which  case  there  are  no  pliable  dements 
that  are  able  to  distort  and  thus  allow  for  a  smooth  variation  of  mesh  displacement. 

Figure  2*3-9  (a,b,c)  shows  the  same  portions  of  the  mesh  for  £=0.8  which  is  an  intermediate  value  of  the 
stiffness  parameter  (0.0<x<2.0).  The  close-ups  of  the  mesh  in  the  vicinity  of  the  leading  and  trailing  edge 
show  that  while  the  mesh  is  not  as  high  a  quality  as  for  it  is  of  sufficient  quality  as  to  allow  an 

accurate  solution.  In  addition,  the  lower  stiffness  parameter  value  has  allowed  larger  grid  elements  between 
the  airfoil  and  the  lower  wall  to  distort  in  a  smooth  fashion,  thus  preventing  the  delamination  phenomenon 
associated  with  the  higher  stiffness  parameter. 
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Figure  2*3-8*  Gdd  for  %=  I +5  when  the  airfoil  is  located  at  its  maximum  vertical  displacement. 
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2.4  Validation  of  Loci-STREAM  with  Experimental  Data  for  Pure  Plunging 
SD  7003  Airfoil 


Experimental  data  (McGowan  et  al.  2008)  involving  phase-averaged  measurements  using  particle  image 
velocimetry  in  a  water  tunnel  was  used  to  validate  the  Loci-STREAM  code  for  an  oscillating  SD7003  airfoil 
in  the  pure-pitch  mode.  Reynolds  numbers  of  10,000  and  40,000  (based  on  free  stream  velocity  and  airfoil 
chord)  were  selected  as  representative  values  for,  respectively,  a  situation  where  transition  would  not  occur, 
and  a  situation  where  transition  in  attached  boundary  layers  would  be  of  some  significance.  The  plunging 
motion  is  at  the  reduced  frequency  of  k=3.93  with  the  pivot  point  at  the  quarter  chord.  The  computational 
and  experimental  results  are  plotted  side-by-side  in  Figure  2.4-1  and  Figure  2.4-2  for  Re=I0,000  and  in 
Figure  2.4-3  and  Figure  2.4-4  for  Re=40,000.  Contours  of  streamwise  velocity  (non-dimensionalized  by  the 
free  stream  velocity)  are  plotted  in  Figure  2.4-1  and  Figure  2.4-3.  Contours  of  vorticity  (non-dimensionalized 
by  free  stream  velocity  and  airfoil  chord)  are  plotted  in  Figure  2.4-3  and  Figure  2.4-4.  The  contours  are 
plotted  at  four  different  phases  at  which  the  data  was  taken,  namely,  the  beginning  of  the  cycle,  %  of  the 
cycle,  Vi  of  the  cycle,  and  3/«  of  the  cycle. 


Exp.  o=3  4 


Loci-STREAM  (SST  model) 


Figure  2.4-1.  SD7003  airfoil;  pure  plunge;  Re= 10,000.  Streamwise  velocity. 
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Loci-STREAM 


Figure  2.4-2.  SD7003  airfoil;  pure  plunge;  Re= 10,000.  Vorticity. 
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Exp.  o=0 


Exp.  Q“1  4 


Exp.  o=3  4 


0.1_ QJ> DJ 


-0.1  0  3  0.6  1.0 


Loci-STREAM 


Figure  2.4-3.  SD7003  airfoil;  pure  plunge;  Re=40,000,  Streamwise  velocity. 
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Figure  4.  SD7003  airfoil;  pure  plunge;  Re=40,00G.  Vorticity. 
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Loci-STREAM 


Figure  2.4-4.  SD7003  airfoil;  pure  plunge;  Re-40,000.  Streamwise  velocity. 


The  agreement  between  the  computation  and  the  experiment  is  very  good.  The  computations  accurately 
capture  the  strong  shear  layer  that  develops  near  the  trailing  edge  during  the  plunge  motion.  The  vorticity 
results  from  the  Loci-STREAM  code  show  its  capacity  to  adequately  resolve  shed  vortical  structures  near  the 
airfoil  and  in  the  wake. 
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2*5  Development  of  an  Optimization  Framework 

In  order  to  create  an  effective  tool  for  the  design  of  MAVs,  the  basic  fluid-structure  solver  must  be  coupled 
with  an  optimization  strategy.  While  the  fluid-structure  solver  can  provide  detailed  information  about  the 
flow  physics  of  a  particular  MAV  configuration  under  any  given  flight  conditions,  a  complete  parametric 
study  involving  numerous  configurations  with  multiple  flight  conditions  must  be  conducted  in  order  to 
determine  the  potential  “best  design”  for  an  MAV.  Given  the  current  state  of  computing  power,  primarily  in 
regards  to  unsteady  simulations  over  complex  3-D  geometries,  one  must  find  a  way  to  determine  the 
sensitivity  of  the  design  to  the  input  parameters  while  conducting  the  minimum  number  of  CFD  simulations. 
An  effective  optimization  strategy  can  reduce,  by  orders  of  magnitude,  the  number  of  simulations  required 
for  the  complete  characterization  of  the  design  with  respect  to  geometric  and  flow  parameters.  Development 
of  an  optimization  framework  was  a  key  element  of  Phase  2  work. 

While  the  ultimate  goal  of  this  effort  is  to  seek  the  3D  kinematics  yielding  the  best  performance  for  flapping 
wings  with  varying  sizes,  frequency,  structural  complexity,  and  capable  of  handling  forward  flight,  hovering 
and  wind  gust  scenarios,  the  task  presented  in  this  section  deals  with  a  better  understanding  of  the  interaction 
and  effects  of  the  unsteady  flow  mechanisms  and  comparing  potential  kinematic  combinations  with  the  use 
of  surrogate  models.  Reynolds  number  effects  have  been  examined  previously,  and  one  such  consequence, 
the  asymmetric  forward  and  back  stroke,  is  seen  as  a  consequence  of  interactions  with  the  jet  like  flow 
feature.  The  sinusoidal  hovering  kinematics  utilized  here  been  used  in  former  studies  as  well,  and  yet  there 
are  still  questions  even  in  the  simplified  2D  domain.  What  constitutes  “good”  kinematics?  This  is  a  context 
specific  question,  e.g.,  the  measure  of  merit  may  be  lift,  but  then  the  next  relevant  question  is  why,  and  at 
what  cost  (power  consumption).  Why  are  certain  combinations  of  variables  better  than  others  and  can  any 
general  trends  be  stated?  This  is  merely  a  starting  point.  After  a  stronger  foundation  has  been  set  other 
aspects  can  be  addressed.  Transitional  effects  are  relevant  at  Reynolds  numbers  seen  by  MAVs  and  are  one 
of  the  open  challenges  in  the  field.  The  3D  effects  at  relatively  low  aspect  ratios  are  quite  important.  The 
traditional  wing  tip  vortices  are  present  but  there  are  also  highly  complex  flow  field  interactions  which  can 
be  seen  in  computations  of  a  flapping  hawk  moth.  All  of  these  are  possible  avenues  that  can  be  eventually 
integrated  into  the  current  framework. 


Surrogate  modeling  is  one  of  the  optimization  methodologies  used  in  engineering  environments.  Queipo  et 
al.  (2005)  present  an  overview  and  highlights  the  strengths  and  issues  in  using  surrogate  based  analysis  while 
Goel  et  al.  (2006)  specifically  address  shortcomings  of  the  experimental  designs.  Surrogate  modeling 
provides  an  efficient  method  for  mining  information  from  limited  data  sets  which  is  usually  expensive,  be  it 
in  computational  or  experimental  costs.  Examples  of  engineering  applications  include  shape  optimization 
using  response  surfaces  as  well  as  other  surrogate  models  in  the  design  of  rocket  injectors  and  supersonic 
turbines  (Shyy  et  al.  2001).  It  is  seen  that  alternatives  to  gradient  based  optimizations  are  needed,  and  these 
examples  provide  empirical  evidence  of  the  utility  in  using  surrogate  models.  In  the  present  study,  they  are 
used  to  get  an  idea  of  how  the  design  space  behaves  away  from  the  known  data  points.  Not  only  does  this 
help  clarify  the  general  trends,  but  when  the  design  of  experiment  (DOE)  is  done  properly,  the  result  is  more 
effic  ient  use  of  computational  resources. 

2.5.1  TheodorseiTs  Solution 

Theodorsen  (1935,  1942),  is  able  to  separate  the  instantaneous  and  delayed  reactions  of  an  incompressible 
fluid  caused  by  an  unsteady  body.  The  assumptions  include  potential  flow,  small  oscillations,  and  thin  airfoil 
simplifications.  Furthermore,  a  Kutta  condition  is  imposed  on  the  trailing  edge,  and  a  planar  wake  is 
assumed.  This  means  that,  strictly  speaking,  this  solution  doesn’t  include  the  effects  of  vortex  roll-up, 

shedding,  as  well  as  open  separations.  For  pitch,  and  plunge  described  as  a(t)  =  aae[s**2^  and 
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h[t)  =  h^e^1  2ir^ h{t)  —  hael^t+2n<$K  respectively,  the  lift  is  predicted  as 


L  =  7ipb2[h+ Ujx -badj+ 2npUJ>C ( k )  h+Ujx+b^~-ct^d 


(2.9) 


where  b  is  semichord  length,  V ^  is  the  freestream  velocity,  k  is  the  reduced  frequency  defined  asit  =  — 
k  —  — and  C(k )  is  the  Theodorsen’s  lift  deficiency  function*  The  TheodorseiTs  lift  deficiency  function  is  a 

”lO 

complex  function  which  can  be  expressed  using  Hankel  functions  of  first  and  second  kind.  Furthermore,  if 
there  is  a  constant  initial  angle  of  attack,  then  this  steady-state  contribution  can  be  superimposed  to  Eq.  (2.9). 


The  first  term  in  Eq.  (2.9)  is  referred  as  the  added  mass  term.  This  term  gives  the  instantaneous  reaction  of 
the  flow  to  the  airfoil  motion,  and  the  resulting  pressure  distribution.  The  second  term,  the  circulatory  term, 
that  contains  the  Theodorsen’s  lift  deficiency  function,  represents  the  influence  of  the  wake  vortices.  Since 
the  magnitude  C(k )  is  always  less  than  unity,  this  circulatory  term  brings  reduction  in  the  lift  amplitude,  and 
also  phase  lag  into  the  flow. 

Two  particularly  useful  approximations  to  the  above  form  are  the  steady -state, 


L  =  2xpUj[h+Ua,a], 

(2.10) 

and  q  uas  i  -stead  y  app  ro  x  i  m  at  i  on  s. 

L  =  nph 1  |  h+Urd-  bad ]  +  InpUJb  h  +  U,  a  +  b  ^  -  a  jd  , 

(2.11) 

For  the  steady -state  case,  we  have  the  steady -state  formula  except  that  the  angle  of  attack  has  been  replaced 
by  the  effective  angle  of  attack.  In  the  quasi-steady  approximation,  C(k)  ->  1  as  k  ->  0,  and  the  impact  of  the 
shed  vortices  is  ignored. 

2*5*2  Surrogate  Modeling 

The  surrogate  modeling  techniques  are  used  to  gain  a  better  understanding  of  the  impact  of  the  design 
variables.  The  various  time  histories  give  a  telling  story  as  to  what  is  happening  for  each  combination,  but 
illustrating  the  overall  impact  is  not  straight  forward,  even  when  limited  to  three  degrees  of  freedom. 
Surrogate  modeling  provides  a  way  to  stand  back  from  the  trees  and  see  the  forest  if  you  will.  The  process  is 
split  into  three  main  parts.  The  first  is  constructing  the  design  of  experiment  (DOE),  or  method  for  choosing 
how  many  and  for  which  points  to  run  full  CFD  simulations.  These  then  provide  the  necessary  data  for  the 
fitting  of  the  surrogate  models  which  can  be  used  to  approximate  a  quantity  of  choice  at  arbitrary  points 
within  the  design  space.  The  third  piece  of  the  process  is  a  sensitivity  analysis  which  can  be  used  to  quantify 
the  importance  of  each  design  variable,  and  in  some  cases  eliminate  them  from  consequent  refinement 
iterations. 

2.5,3  Design  of  Experiment  (DOE) 

The  DOE  used  a  face  centered  cubic  design  (FCCD)  (Shyy  et  al.  2001)  and  then  Latin  hypercube  sampling 
(LHS)  (Queipo  et  al.  2005)  to  appropriately  fill  in  the  remainder  of  the  design  space.  A  2nd  order  polynomial 
response  surface  construction  has  (N+l)(N+2)/2  coefficients,  N  being  the  number  of  variables,  and  in 
general,  one  wants  twice  this  many  data  points  for  a  proper  curve  fit.  A  FCCD  design  provides  2n+2N+1 
points:  2X  corner  points,  2N  face  points,  and  one  center  point.  Thus  for  three  design  variables,  FCCD 
provides  15  of  the  20  points  required.  The  LHS  then  provides  a  method  for  efficiently  choosing  the  rest  of 
the  points  by  maximizing  the  distance  between  the  added  points,  though  by  no  means  is  it  the  only  alternative 
(Queipo  et  al.  2005). 
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2.5.4  Composite  Surrogates 

Which  surrogate  model(s)  to  use  is  an  interesting  challenge  in  and  of  itself.  The  current  approach  uti  lizes  the 
current  state  of  the  art  work  of  Viana  et  al.  (2008),  Here  a  number  of  different  surrogate  models  are 
constructed  and  are  evaluated  based  on  their  respective  cross-validation  errors,  namely  PRESS,  in  this 
fashion  no  extra  test  points  are  needed,  rather  the  models  are  constructed  with  one  less  training  point  and  the 
deviation  of  this  point  from  the  constructed  surrogate  is  used  calculate  one  component  of  the  PRESS.  Each 
of  the  training  points  is  treated  in  a  similar  manor  and  one  can  subsequently  quantify  how  well  the  respective 
surrogate  model  fits.  As  Viana  et  al,  (2008)  illustrates  for  problems  in  lower  dimensions,  using  the  best 
PRESS  surrogate  model  might  be  justified,  whereas  in  higher  dimensions  it  is  much  riskier  to  do  so. 


The  current  setup  evaluates  Kriging  (Queipo  et  al.  2005),  2nd  order  polynomial  response  (PRS)  (Myers  and 
Montgomery  2002),  radial  basis  neural  network  (RBNN)  (Cheng  and  Titterington  1994),  and  6  support 
vector  regression  (SVR)  (Smola  and  Scholkopf  2004)  models  noted  in  Table  2.5-1  below.  The  two  best 
models  in  the  current  context,  as  measured  by  those  exhibiting  the  lowest  PRESS  values,  are  the  SVR  model 
using  a  full  spline  kernel  and  Kriging. 


Table  2,5-1,  Defining  surrogate  model  traits. 


Model 

Comment  1 

Comment  2 

Kriging 

Linear  Regression  Model 

Gaussian  Correlation  Model 

PRS 

2nd  Order  Polynomial 

RBNN 

Max  Neurons  =  1000 

SVR  1 

Linear  Spline  Kernel 

Full-  infinity  as  upper  bound 
(non-separable  case) 

SVR2 

Linear  Spline  Kernel 

Short-  finite  upper  bound 

SVR3 

Exponential  Kernel 

Full 

SVR4 

Exponential  Kernel 

Short 

SVR5 

Gaussian  Kernel 

Full 

SVR6 

Gaussian  Kernel 

Short 

2.5.5  Global  Sensitivity  Analysis 

The  global  sensitivity  analysis  is  in  general  useful  for:  (i)  Determining  if  a  variable  is  particularly  influential 
in  the  design  space,  if  not  perhaps  the  variable  can  be  fixed  and  the  degrees  of  freedom  and  complexity  of  the 
problem  reduced,  (ii)  Ranking  the  importance  of  the  design  variables,  (iii)  Quantifying  the  degree  of 
coupling  between  design  variables.  For  example,  is  the  influence  on  the  design  space  mostly  an  individual 
effort,  or  is  there  an  effect  caused  by  the  interaction  of  variables? 

Sobol’s  method  (Sobol  1993)  is  used  to  for  the  global  sensitivity  evaluations.  The  surrogate  model  can  be 
written  as: 

/(*)  =  /» +£/(*,)+ £//(*<  >*,)+•••>  (2-12) 

* 

Once  this  decomposition  has  been  calculated  the  total  variance, 

D=  \f{x)dx-f* ,  (2.13) 
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and  partial  variances,  e.g.? 

A  =  \fi{x)dx\dxidx^  (2.14) 

can  be  calculated.  In  this  fashion,  individual  contributions,  such  as  D|/D,  or  combinations  of  variables,  e,g., 
D|2/D,  can  be  quantified,  effectively  capturing  the  sensitivity  of  the  variable(s)  under  consideration. 


2,5,6  Elliptical  Airfoil 

A  plunging  and  pitching  airfoil  was  used  for  this  purpose.  The  translational  and  rotational  airfoil  motions  are 
dictated  by  Eqs.  (2. 1 5)and  (2. 1 6): 

h(t)  =  has\n{2xft),  (2.15) 

a(t)  -a0 -aa  sin(2xfl  +  0) ,  (2.16) 

-o  h 

- -  2ti. - ■ 

Figure  2.5-1.  Illustration  of  the  kinematic  parameters  for  normal  hovering. 

Here  h(t)  and  ha  are  the  translational  position  and  plunging  amplitude  respectively  (see  Figure  2,5-1).  The 
angular  orientation,  initial  angle,  and  angular  amplitude  are  a(t ),  a q*  and  a£,.  The  pitching  is  about  the  center 
of  the  rigid  airfoil;  this  is  an  el  lipse  having  a  15%  thickness  for  all  cases  under  consideration.  The  phase  lag 
between  the  two  motions  is  t|>,  and  the  frequency  is  denoted  /whereas  the  time  is  again  /.  While  there  are  a 
few  choices  in  how  to  accommodate  these  kinematics  computationally,  the  current  implementation  forces  the 
grid  to  rotate  and  translate  with  airfoil.  The  geometric  conservation  law  (GCL),  a  necessary  consideration  in 
domains  with  moving  boundaries,  is  satisfied.  The  boundary  condition  applied  to  all  outer  boundaries  is  the 
incompressible  inlet  with  density  and  velocity  specified. 


Due  to  the  kinematic  constraints  there  are  only  two  relevant  non-dimensional  groups  in  the  incompressible 
case.  The  Reynolds  number  is  given  by 


Re-- 


nf^ref  (2^A)£ 


U„,L, 


V  v 

and  the  plunging  amplitude  to  chord  ratio  is  given  by 

2h 


(2.17) 


(2.18) 


The  reduced  frequency  is  given  by 


^Lni  2jijc 

2Unf  2(2*A) 


(2.19) 


One  can  see  that  the  reduced  frequency,  k  expressed  in  Eq.  (2.19),  contains  the  same  information  as  the 
plunging  ratio.  Note  that  since  the  Reynolds  number  is  held  constant,  where  the  reference  velocity  and 
length  are  taken  to  be  the  maximum  translational  velocity  and  the  chord  length  respectively,  ha  and  f  are  not 
independent.  The  three  quantities  that  can  vary  are  hm  o*j,  and  tp.  For  the  compressible  case  the  Mach  number 
is  also  relevant. 
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2.5.6.1  Design  Space 

The  range  of  the  variables  was  chosen  based  on  the  works  of  Weis-Fogh  (1973)  and  Ellington  (1999),  which 
tabulated  actual  measurements  from  various  species,  some  of  which  are  reproduced  in  Table  2.5-2.  The 
bounds  for  the  simulations  are  listed  in  Table  2.5-3,  though  it  is  worth  noting  that  (p  <  90*,  also  referred  to  as 
delayed  rotation,  is  not  expected  to  be  found  in  nature  for  hovering  flight  as  the  airfoil  is  flying  upside  for 
portions  of  the  cycle*  These  cases  were  included  to  provide  some  symmetry  for  the  other  extreme,  cp  >  90°  or 
advanced  rotation.  Furthermore,  measurements  regarding  aa  are  also  difficult  to  come  by,  though  comments 
found  in  Ellington  (1999)  would  imply  that  the  choice  of  bounds,  Table  2*5-3,  for  this  variable  are 
reasonable. 


c  (cm) 

f  (Hz) 

2hjc 

Pawing 

Fruit  Fly: 
Drosophila  virilis 

0.15 

240 

3.5 

250 

Honey  Bee: 

Apis  mellifica 

0.43 

240 

2.8 

1900 

Bumble  Bee: 

Bombus  terrestris 

0.73 

156 

2.8 

4800 

Hummingbird: 

Archilochus 

colubris 

1.5 

52 

3.6 

6400 

Hawkmoth: 

Manduca  Sexta 

2.5 

27.3 

2.6 

6700 

Hummingbird: 

4.3 

15 

3.6 

15000 

Patagona  gsgas 

Table  2.5-2.  Selected  data  on  the  time  and  length  scales  encountered  in  nature.  The  examples  listed  do  not 
provide  upper  or  lower  bounds  for  any  of  the  categories  listed,  but  do  provide  a  window  which  captures  many  of 
the  insects  and  animals  capable  of  hovering  flight. 


Parameter 

Minimum 

Maximum 

2Vc 

2.0 

4.0 

a. 

45’ 

80* 

<D 

60’ 

120* 

Table  2.5-3.  Minimum  and  maximum  values  of  the  plunging  ratio,  angular  amplitude,  and  phase  lag  simulated 
in  the  surrogate  modeling  exercise. 


To  balance  the  computational  expense  these  simulations  were  carried  out  on  an  81x81  grid,  a  resolution 
which  was  not  grid  independent,  but  captured  the  relevant  behavior  immediately  surrounding  the  airfoil  fairly 
well. 

2.5.6*2  Force  Interpretation 

To  better  understand  the  implications  and  [imitations  of  the  surrogate  modeling  results  an  example  is 
presented  of  a  representative  normal  hovering  case  at  a  Re  of  100  (Figure  2*5-2)*  The  objective  functions  for 
these  cases  will  be  the  time  integrated  CL  and  an  approximation  to  the  non-dimensional  time  averaged  power 
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required  which  is  calculated  by  multiplying  the  non-dimensional  translational  velocity  by  the  CD. 


The  discussion  following  is  generally  applicable  to  the  cases  where  0^90,  with  slight  modifications.  For 
cases  where  the  phase  lag  dictates  advanced  rotation  (0  >  90)  or  delayed  rotation  (<X>  <  90),  the  same  ideas 
can  be  extended  though,  like  the  parameter  suggests,  the  translation  and  rotation  will  be  out  of  phase.  The 
cycle  can  be  broken  up  info  three  overlapping  regions  defined  by  the  unsteady  flow  mechanisms  present. 
The  first  region  starts  at  point  1,  once  again  referring  to  Figure  2.5-2,  which  is  near  a  local  minimum  in  the 
lift.  As  the  airfoil  is  vertical  at  This  point,  one  would  generally  expect  zero  lift.  As  time  continues  the  airfoil 
turns  back  into  its  previous  trajectory  which  is  commonly  referred  to  as  wake  capturing,  points  1,  2,  and  3, 
The  peak  seen  at  point  3  will  be  referred  to  as  the  wake  capturing  peak.  Flow  field  snapshots  of  vorticity 
(Figure  2.5-3),  demonstrate  the  nomenclature  more  clearly. 


Figure  2.5-2.  Illustration  of  the  force  histories  over  a  normal  hovering  cycle  for  the  parameters  of  2hB/c  —  3,0, 
a*  =  45,  and  0=90  and  the  corresponding  airfoil  positions. 
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Figure  2,5-3*  Illustration  of  vorticity  couutours  (red:  clockwise,  blue  counter-clockwise)  and  the  wake 
capturing  portion  of  the  cycle  where  the  airfoil  turns  back  onto  its  wake  possibly  taking  advantage  of  the 
unsteady  aerodynamics.  In  Figure  2.5-2,  these  plots  would  correspond  to  points  on  either  side  of  point  1.  The 
unlabelled  point  0  on  the  left  which  belongs  to  the  end  of  the  backstroke,  and  point  2  on  the  right* 


The  second  and  third  unsteady  flow  features  overlap  significantly*  The  most  commonly  known  is  the  delayed 
stall  phenomenon  which  is  generally  associated  with  a  leading  edge  vortex  (Figure  2.5-4).  Here  a  vortex 
forms  behind  the  leading  edge  of  the  airfoil,  causing  a  low  pressure  region  and  enhancing  lift.  Note  that  in 
the  case  illustrated,  higher  l  ift  is  achieved  at  angles  of  attack  of  45  degrees,  an  angle  well  beyond  the  steady 
state  stalk  In  cases  with  higher  angular  amplitudes,  and  therefore  lower  angles  of  attack,  the  peak  at  points  1 
and  8  can  be  reduced  significantly  because  the  orientation  of  the  airfoil  is  not  able  to  promote  LEV 
formation. 


Figure  2.5-4.  Illustration  of  unsteady  delayed  stall  mechanism  and  the  LEV  with  vorticity  countours  (red: 
clockwise,  blue:  counter-clockwise)  In  Figure  2*5-2,  these  snapshots  correspond  to  points  4,  6,  and  8 
respectively. 
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Figure  2*5-5*  Illustration  of  the  jet  flow  feature  which  develops  in  the  path  of  the  airfoil.  Velocity  contours  are 
displayed  (red:  positive,  blue:  negative)  with  select  percentages  of  the  maximum  translational  velocity 
marked*  In  Figure  4  these  snapshots  would  correspond  to  points  2 , 6,  and  10  respectively* 


The  last  unsteady  flow  mechanism  is  a  jet- 1  ike  flow  feature  (see  Figure  2,5-5)  present  in  the  path  of  the 
airfoil  where  a  persistent  downward  velocity  develops.  These  regions  within  the  jet  are  of  comparable 
magnitude  to  the  maximum  translational  velocity  of  the  airfoil  itself  and  are  influential  for  a  large  segment  of 
the  cycle,  in  this  case  roughly  from  point  4  to  point  9.  This  jet  is  the  cause  of  the  asymmetry  seen  between 
forward  and  backstrokes  (see  Figure  2.5-3),  and  is  due  to  the  fact  that  the  jet  is  slightly  off  centered  in  the 
direction  in  which  the  motion  was  started.  The  jet  also  explains  the  local  minimum,  near  point  5.  If  the  angle 
of  attack  is  low  (higher  angular  amplitudes),  or  the  jet  strength  is  stronger  (shorter  plunging  amplitudes  mean 
the  jet  decays  less  between  encounters)  then  this  wake  valley  will  become  deeper.  To  help  convince  the 
reader  of  the  jet’s  influence  CL  has  been  plotted  during  the  first  two  cycles  as  well  as  the  1 5lh  cycle  where  the 
differences  between  previous  cycles  has  largely  stopped  at  the  spatial  and  temporal  resolutions  used  (Figure 
2.5-3).  Note  that  because  the  simulation  starts  with  a  targe  discontinuity  in  airfoil  velocity,  and  the  plots  are 
shifted  such  that  the  force  histories  start  when  the  airfoil  is  at  the  end  of  its  translation,  the  1st  cycle  does  not 
imply  no  wake.  Rather  it  serves  to  provide  credibility  to  the  claim  that  the  jet  is  at  least  partly  responsible  for 
the  aerodynamic  performance  recorded  (i,e.  the  LEV  is  not  the  only  factor)  as  it  is  clear  that  as  the  jet 
changes  strength,  the  aerodynamic  characteristics  respond  in  a  noticeable  fashion. 
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Figure  2,5-6.  The  lift  coefficient  (top)  and  power  required  (bottom)  as  approximated  by  the  WAS  illustrating 
the  effect  of  the  plunging  amplitude.  For  the  integrated  effect  over  the  entire  cycle  it  is  seen  that  this  is  the 
least  sensitive  of  the  design  variables.  Arrows  show  the  direction  of  increasing  plunging  amplitude. 


2.5.6.3  Kinematics 

Plunging  Amplitude 

Somewhat  surprisingly,  the  plunging  amplitude,  and  with  it  the  reduced  frequency,  had  the  smallest  influence 
on  lift  and  the  required  power  throughout  the  study.  Since  the  peak  translational  velocity  of  the  airfoil  is 
equal  for  all  cases,  the  fact  that  the  plunging  amplitude  has  little  effect  on  the  power  required  per  cycle  is  to 
be  expected.  However,  the  effect  on  lift  of  the  unsteady  flow  mechanisms  is  less  easily  predicted  as  is  their 
cumulative  impact.  It  is  readily  apparent  from  the  surrogate  response  of  the  design  space,  as  seen  by  the 
small  gradients  found  in  the  direction  of  the  arrows  in  Figure  2.5-6,  that  the  plunging  amplitude’s  influence  is 
limited. 
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To  gain  a  more  comprehensive  look  on  the  plunging  amplitude’s  influence,  three  force  histories  are  shown  in 
Figure  2.5-7.  The  first  observation  is  that  the  wake  capturing  peak  is  slightly  smaller  as  the  plunging 
amplitude  is  increased.  The  most  significant  change  though  is  in  the  wake  valley.  The  valley  is  the  deepest 
when  the  plunging  amplitude  is  short.  As  the  plunging  amplitude  increases,  the  intensity  of  the  jet  found  in 
the  wake  is  able  to  decay  before  interaction  with  the  airfoil.  The  influence  of  the  jet  is  detrimental  to  the  lift 
under  these  conditions,  and  a  weaker  jet  corresponds  to  a  higher  lift  for  airfoils  at  identical  angles  of  attack 
and  translational  velocities.  Moving  on  to  the  delayed  stall  peak,  seen  as  the  global  maxi  mums,  there  is  little 
influence  on  the  peak  amplitude.  The  general  implication  then  is  to  minimize  the  negative  impact  the 
persistent  jet  has  after  the  wake  capturing  peak.  In  the  context  of  plunging  amplitude  this  means  increasing 
the  stroke  length  to  encounter  a  weaker  jet,  or  looking  at  it  from  another  perspective  this  could  be  seen  as 
approaching  jets  of  similar  intensity  at  higher  angles  of  attack. 


Figure  2.5-7.  Snapshots  of  the  effect  of  the  plunging  amplitude  (2.0,  3.0,  4.0)  on  the  instantaneous  force 
history. 

Angular  Amplitude 

In  direct  contrast  to  the  discussion  concerning  the  plunging  amplitude,  the  design  space  is  quite  sensitive  to 
the  angular  amplitude,  see  Figure  2.5-8.  In  general  it  is  found  that  lower  angular  amplitudes,  thus  higher 
angles  of  attack,  lead  to  higher  power  requirements  and  lift.  The  power  required  result  is  expected  as  higher 
angles  of  attack  will  correspond  to  more  drag,  and  consequently  a  larger  power  requirement.  The  lift  result, 
while  not  unexpected,  is  the  result  of  interactions  of  the  unsteady  aerodynamics  and  is  not  entirely  intuitive. 
What  is  evident  is  the  fact  that  higher  lift  and  lower  power  required  aspirations  are  in  direct  competition  on 
either  end  of  the  angular  amplitude  spectrum. 
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Figure  2.5-8.  The  lift  coefficient  (top)  and  power  required  (bottom)  as  approximated  by  the  WAS  illustrating 
the  effect  of  the  angular  amplitude.  A  great  deal  of  variance  occurs  with  these  changes  and  it  is  seen  that  in 
general  lift  coefficient  over  the  whole  cycle  benefits  from  lower  angular  amplitudes.  Arrows  show  the 
direction  of  increasing  angular  amplitude. 

Looking  again  at  force  histories,  the  focus  this  time  is  angular  amplitude,  see  Figure  2.5-9.  The  first 
observation  is  that  the  amplitude  of  the  wake  capturing  peak  does  not  change  much  with  angular  amplitude. 
This  is  somewhat  surprising  considering  that  higher  angular  amplitudes  lead  to  higher  angular  velocities 
which  then  interact  with  the  wake  behind  the  airfoil  at  the  end  of  the  translation.  The  root  of  why  the  angular 
amplitude  matters  so  much  is  seen  clearly  in  the  wake  valley  and  subsequent  delayed  stall  peak.  As  the 
angular  amplitude  is  increased,  and  the  angles  of  attack  decreased,  the  airfoil  approaches  the  jet  at  less 
favorable  orientations,  making  the  wake  valley  deeper  and  wider.  However,  the  delayed  stall  peak  also 
decreases  and  so  we  have  two  effects  acting  in  concert  The  lower  angles  of  attack  experienced  do  not 
provide  diminished  (and  sometimes  no)  support  for  LEV  formation  and  thus  a  much  lower  maximum  lift 
value.  At  the  same  time,  the  higher  angular  amplitudes  also  provide  a  more  streamlined  body,  producing  less 
drag.  It  is  seen  that  the  rule  of  thumb  is  higher  angles  of  attack  and  lower  angular  amplitudes  provide  higher 
lift  through  a  combination  of  jet  interaction  and  LEV,  but  at  the  same  time  create  a  higher  drag. 
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Figure  2.5-9.  Snapshots  of  the  effect  of  angular  amplitude  (45*  62. 5*  80)  on  the  instantaneous  force  histones. 
Phase  Lag 


For  the  final  design  variable*  phase  lag,  the  general  trend  is  for  O  >  90,  or  advanced  rotation,  to  produce 
higher  lift  values  (Figure  2.5-10).  How  this  occurs  will  be  explained  shortly  but  is  once  again  a  combination 
of  unsteady  effects.  It  is  also  interesting  to  note  that  the  phase  lag's  influence  on  the  integrated  lift  is  minimal 
when  angular  amplitudes  are  high.  Yet,  at  higher  angular  amplitudes  is  exactly  where  largest  influence  is  felt 
when  considering  the  power  required.  The  fact  that  the  most  efficient  mode,  from  an  energy  input 
perspective*  is  when  <!>  «  90.  This  puts  the  minimum  angle  of  attack  at  the  maximum  translational  velocity* 
and  therefore  less  drag  and  power  required. 


Coupled  with  the  force  histories  in  Figure  2.5-11*  a  more  complete  picture  emerges.  During  <D  <  90,  or 
delayed  rotation,  the  lift  starts  off  negative  and  the  wake  capturing  peak  is  shifted.  From  the  illustrations  of 
the  airfoil  positions*  it  can  be  seen  that  the  airfoil  is  temporarily  flying  upside  down  and  hence  the  negative 
lift  occurs  immediately  upon  changing  direction.  When  the  wake  capture  peak  does  occur  it  is  also  being 
influenced  by  the  jet,  and  at  the  higher  angles  of  attack  the  interaction  is  beneficial.  The  delayed  stall  peak 
suffers.  Though  there  is  a  higher  angle  of  attack  occurring  at  the  maximum  translational  velocity*  as 
compared  to  O  =  90,  the  angular  velocity  is  negative  and  not  conducive  to  LEV  formation. 
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Figure  2.5-10.  The  lift  coefficient  (top)  and  power  required  (bottom)  as  approximated  by  the  WAS  illustrating 
the  effect  of  the  phase  lag.  Again  a  pronounced  response  in  most  areas  of  the  design  space.  In  general  higher 
phase  lags,  corresponding  to  advanced  rotation,  are  more  beneficial  from  a  lift  point  of  view  though  that 
statement  is  not  universally  true,  especially  at  higher  angular  amplitudes.  Arrows  show  direction  of 
increasing  phase  lag. 


For  0  >  90,  or  advanced  rotation,  the  airfoil  starts  turning  earlier  such  that  a  positive  angle  of  attack  is 
achieved  upon  wake  capture,  producing  favorable  lift.  This  is  immediately  followed  by  a  pronounced  wake 
valley.  It  was  seen  that  when  <E>  =  90,  the  wake  valley  could  severely  impede  lift.  That  effect  is  exaggerated 
here  where  lower  angles  of  attack  encounter  the  jet.  After  that  interaction,  a  very  favorable  delayed  stall  peak 
occurs  as  a  higher  angle  of  attack,  with  positive  angular  velocity,  is  present  at  the  maximum  translational 
velocity.  While  one  can  see  from  the  surrogates  the  consequences  of  changing  the  phase  lag,  the  variable's 
influence  is  more  than  initially  implied.  As  seen  from  the  force  histories,  the  influences  on  the  wake  capture, 
jet  interaction,  and  delayed  stall  partially  cancel  out  thus  obscuring  its  overall  importance. 
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Figure  2.5-11.  Force  histories  and  airfoil  positions  illustrating  delayed  (left)  and  advanced  (right)  rotation. 
Here  we  see  the  utility'  in  the  advanced  rotation  as  it  keeps  the  lift  levels  high  after  the  influence  of  the  delayed 
stall  has  subsided. 


2.5.64  Sensitivity  Analysis 


Figure  2.5-12.  Sensitivity  indices  for  <C(>  (left)  and  <Prttl>  (right)  in  normal  hovering  (no  freestream).  The 
main  and  total  variances  are  given  from  left  to  right:  plunging  amplitude,  angular  amplitude,  phase  lag. 


There  is  some  sense  ofimportance  conveyed  by  the  plots  of  the  design  space  in  Figure  2.5-6,  Figure  2.5-8 
and  Figure  2.5-10.  In  low  dimensions  it  may  not  be  difficult  to  discern  relative  importance  amongst  the 
design  variables.  However,  Figure  2.5-12  provides  a  more  objective  measure  of  influence.  The  global 
sensitivity  indices  are  tabulated  and  a  more  quantifiable  relationship  can  be  determined.  As  expressed  earlier, 
depending  on  whether  instantaneous  effects  are  important  or  whether  the  integrated  result  is  sufficient,  the 
role  of  phase  lag  could  be  underestimated.  Regardless,  within  the  design  space  examined,  the  plunging 
amplitude  (and  therefore  the  reduced  frequency)  plays  a  small  role  in  determining  the  unsteady  flow  physics 
seen  by  an  airfoil  in  normal  hovering. 


35 


STTR  Phase  2:  A  Simulation  Environment  for  Aerodynamic  Analysis  &  Design  of  Flapping  Wing  Micro  Air  Vchides-Final  Report 


Figure  2.5-13,  Vorticity  and  vertical  velocity  contour  plots  of  normal  hovering  at  Re=100  (left)  and  Re=1000 
(right) 


2.S.6.5  Reynolds  Number 

As  the  Reynolds  number  is  increased,  the  vortex  dynamics  quickly  lead  to  chaotic  behavior.  As  seen  in 
Figure  2.5-13,  stronger  vortices  are  created  as  the  Reynolds  number  increases.  These  vortices  persist  and 
interact  with  each  other  over  multiple  stroke  cycles  which  in  turn  leads  to  the  chaotic  force  histories 
experienced  by  the  airfoil  (see  Figure  2.5-14).  Another  implication  of  raising  the  Reynolds  number  is  the 
strengthening  of  the  jet-like  flow  feature  which  loses  a  stable  configuration. 


Figure  2.5-14.  Reynolds  number  effect  on  CL  going  from  Re=100  (Red)  which  is  displays  a  repeatable 
pattern,  to  Re=lO00  (Dashed-Blue)  whose  force  histories  vary  substantially  from  stroke  to  stroke* 
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To  summarize,  the  work  described  in  this  section  has  addressed  the  following: 

( 1 )  Gain  a  better  understanding  of  the  2D  flow  physics  present  in  normal  hovering: 

The  influence  of  kinematics,  and  the  Reynolds  number  was  investigated.  Regarding  the  kinematics,  it  was 
found  that  under  most  of  the  circumstances  studied,  the  plunging  amplitude/  reduced  frequency  played  a 
relatively  minute  role  in  determining  the  flow  physics  and  consequent  forces  experienced  by  the  airfoil*  As 
the  reduced  frequency  is  one  of  the  non-dimensional  groups  it  was  somewhat  expected  that  the  parameter 
would  have  a  larger  influence*  A  point  of  caution  is  warranted  as  the  conclusions  presented  are  taken  in  the 
context  of  the  current  design  space*  The  flow  quickly  tended  towards  chaos  due  to  the  increased  intensity  and 
number  of  vortex  interactions  as  the  Reynolds  number  was  increased. 

(2)  Demonstrate  the  utility  of  surrogate  modeling  in  building  a  quantitative  framework  to  assess  hovering 

kinematics  and  the  associated  flow  physics: 

The  WAS  provided  insight  regarding  general  trends  not  easily  extracted  otherwise  (e.g.  examining  the  force 
histories  at  the  DOE  points)*  The  global  sensitivity  evaluations  quickly  show  the  relative  influence  each 
design  variable  carries,  and  one  can  quickly  evaluate  a  variable’s  importance.  In  the  cases  investigated 
above,  the  angular  amplitude  was  always  sensitive  to  small  changes*  As  was  explained,  the  influence  of  the 
phase  lag  was  partially  hidden  by  the  competition  in  effects  in  the  zero  freestream  cases,  so  one  must 
exercise  caution  when  making  more  general  implications  (e.g*  In  the  present  context  of  integrated  values 
over  the  course  of  a  cycle,  the  surrogates  give  what  was  asked  of  them.  If  one  wanted  to  make  conclusions 
about  the  instantaneous  effects,  one  would  have  to  look  more  closely  at  the  flow  physics,  as  was  done,  or 
redefine  a  proper  objective  function).  It  is  also  worth  noting  that  most  of  the  results  show  the  design  variables 
are  largely  uncoupled.  An  implication  of  this  is  that  a  simple  liner  superposition  of  design  variable  effects 
may  suffice  in  the  present  context  In  the  work  conducted  so  far,  it  is  found  that  the  WAS  provides  an 
efficient  method  for  characterizing  the  flow  physics  and  quantifying  kinematic  effects  in  normal  hovering. 


2. 6  Coupling  of  Loci-STREAM  and  IS  LA  MS  for  Fluid-Structure  Interaction 

During  the  last  phase  of  the  project,  work  was  undertaken  to  couple  the  fluid-dynamics  solver  Loci- 
STREAM  developed  at  Streamline  Numerics  Inc.  (SNI)  with  the  structural  solver  NLAMS  (Nonlinear 
Membrane  Shell  Solver)  developed  at  the  University  of  Michigan  (UM).  As  a  first  step,  the  requirements 
necessary  to  integrate  the  two  codes  were  formalized.  The  following  are  some  of  the  key  outcomes  of  the 
collaboration: 

I)  Due  to  the  relatively  small  computational  requirements  of  the  structural  solver  compared  to  the  fluids 
solver,  it  was  decided  that  no  attempt  would  be  made  to  parallelize  the  existing  structural  solver.  Thus,  the 
existing  implementation  of  structural  solver  will  compute  the  entire  structural  dynamics  problem  for  the  wing 
on  each  process,  while  the  fluid  dynamics  solver  will  partition  its  part  of  the  combined  application  over 
multiple  processes,  each  solving  only  its  portion  of  the  domain.  The  only  limitation  of  this  approach  is  the 
requirement  that  the  complete  structural  dynamics  problem  (for  the  whole  wing  for  example)  be  capable  of 
residing  in  memory  on  each  process  without  causing  paging  from  the  system  disk  to  system  RAM.  Due  to  the 
small  size  of  the  structural  dynamics  meshes  currently  being  used  for  the  flapping- wing  problems  (less  than 
10,000  cells),  this  requirement  is  easily  met*  When  larger  meshes  are  used  for  the  structural  dynamics  part  of 
the  calculation  (estimated  to  be  greater  than  several  hundred  thousand  cells),  work  will  have  to  be  undertaken 
to  parallelize  the  structural  dynamics  solver*  In  addition,  given  the  relatively  small  nature  of  the  structural 
dynamics  mesh  (approximately  100  times  smaller  than  the  fluid  dynamics  mesh),  no  detrimental  effect  on  the 
overall  scalability  of  Loci-STREAM  for  the  combined  problem  is  anticipated. 
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2)  An  interface  has  previously  been  established  between  the  structured  solver  STREAM  and  NLAMS.  Since 
both  of  these  solvers  are  written  in  Fortran  90,  this  interface  relied  heavily  on  the  use  of  Fortran  modules  to 
communicate  data  bewtween  the  fluids  and  structural  solvers.  Since  Loci -STREAM  is  written  in  C++  and 
thus  cannot  access  variables  via  modules  from  Fortran  90,  an  interface  wrapper  was  written  for  NLAMS 
which  allows  all  global  information  required  by  NLAMS  to  be  passed  via  function  arguments. 


3)  The  existing  NLAMS  interface  was  restructured  and  substantially  simplified  for  use  with  Loci-STREAM. 
Due  to  the  nature  of  Loci-STREAM  it  was  deemed  difficult  to  implement  the  fluid/structure  iteration  loop 
which  was  used  to  coordinate  the  combined  fluid  and  structures  calculations  between  the  structured  solver 
STREAM  and  NLAMS,  Thus,  it  was  decided  to  dispense  with  the  fluid/structure  loop  and  instead  use  the 
existing  iteration  loop  within  Loci-STREAM  along  with  a  modulo  input  flag  (NLAMS  called  once  for  every 
10  fluid  iterations,  for  example)  to  control  the  coupling  between  the  fluid  dynamics  and  structural  dynamics 
problems. 

With  a  clearly  defined  interface  now  established  between  NLAMS  and  Loci-STREAM,  a  CSD 
(Computational  Solid  Dynamics)  module  is  being  developed  for  Loci-STREAM  which  will  gather  all 
required  data  from  the  fluid  dynamics  mesh  and  solution  required  by  NLAMS  to  solve  the  structural 
dynamics  problem.  In  addition,  a  general  restructuring  of  the  existing  Loci-STREAM  grid  motion  module  to 
support  solution-dependent  mesh  deformation  problems  (problems  in  which  the  mesh  deformation  is  a 
function  of  the  flow  solution  as  opposed  to  simply  a  time-dependent  forced  boundary)  has  recently  been 
completed. 

2.6. 1  Test  Case 

As  preparation  for  testing  and  validation  of  the  fluid-structure  capability  using  Loci-STREAM, 
computational  models  were  set  up  to  simulate  an  aluminium  monolithic  Zimmerman  wing  flapping  at  10  Hz 
with  +/-  17  or  21  degree  amplitude  by  coupling  the  U  M/NLA  MS,  MSCMarc,  and  STREAM  codes.  A  new 
grid  topology  was  developed  for  this  wing  which  helped  overcome  some  of  the  re-meshing  problems  of 
STREAM.  Two  different  views  of  the  CFD  mesh  developed  in  Pointwise  are  given  in  Figure  2,6- L  It  may 
be  observed  that  one  of  the  block  interfaces  previously  at  the  “tip”  is  now  moved  more  inboard.  This  had  a 
major  impact  on  the  negative  Jacobian  problem  which  generally  occurred  in  the  tip  region. 


Computational  meshes  for  fluid  and  structural  simulation  are  built  as  shown  in  Figure  2.6-1  and  Figure  2.6-2, 
respectively.  Details  of  mesh  resolutions  are  shown  in  Table  2.6-1  and  Table  2.6-2.  Note  that  in  order  to 
avoid  problems  in  the  re-meshing  phase  during  the  fluid  computation,  a  slightly  large  first  grid  space  is 
utilized  for  the  '"base”  grid  of  the  CFD  model.  For  the  CFD  mode!,  the  distance  of  the  outer  boundary  from 
the  wing  is  25c  (0.625  mm)  and  the  thickness  of  wing  is  zero.  For  the  CSD  model,  the  model  called  'base' 
does  not  include  a  square  rigid  region  near  the  wing  root,  whereas  the  models  named  'medium1  and  ‘fine1 
include  this  region.  The  thickness  of  the  structure  is  0,4  mm,  and  Young's  modulus  and  density  of  the 
material  (aluminium)  are  70.0  GPa  and  2700  kg/m",  respectively.  Two  wing  kinematics  are  considered:  (a) 
1  —  cos  (a >£);  (b)  sin(a)t),  where  a>  is  the  angular  velocity  of  flapping  motion.  Flapping  frequencies  are  10 
and  15  Hz.  Flapping  amplitudes  are  17  and  21  deg.  All  parameters  and  dimensionless  parameters  are 
summarized  in  Table  2.6-3. 
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Plane  Overall 

Figure  2,6-1.  Images  of  computational  mesh  for  fluid  dynamics. 


Medium  Fine 

Figure  2*6-2.  Images  of  structural  meshes. 


Table  2.6-1.  Computational  meshes  for  fluid  dynamics. 


Total  number  of  cells 

First  grid  spacing 

Base 

0.7  Million 

2.5e-3  (0.1  c) 

Medium 

0.7  Million 

l.0e-3  (0.04  c) 

Fine 

1.2  Million 

5.0e-4  (0.02  c) 

*  c  indicates  length  of  wing  root  (0.025  mm) 

Tabic  2.6-2.  Computational  meshes  for  structural  dynamics 


Total  number  of  elements 

Root 

Base 

244 

cantilevered 

Medium 

480 

considered  fixed  root 

Fine 

774 

considered  fixed  root 
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Table  2.6-3.  Geometric*  kinematics*  and  material  parameters  of  UF's  Zimmerman  wing 


Name 

Nomenclature 

Number  &  Unit 

Wing  length 

b 

75  mm 

Chord  length  of  wing  root 

Grooi 

25  mm 

Aspect  ratio 

AR-b2/S 

7.65 

Flapping  frequency  (input) 

f 

10-15  Hz 

Flapping  amplitude  (input) 

±17°  (34°)  or  ±21°  (42°) 

±0.6 1 08  rad  ( 1 .22 1 6  rad) 

Young’s  modulus  of  material 

f^alum 

70.0  GPa 

Density  of  material 

Tallinn 

2700  kg/m3 

Reference  velocity 

Uref 

0.89-1 .65  m/s 

Reynolds  number 

Re 

1483-2750 

Reduced  frequency 

k 

0.714-0.882 

Density  ratio 

PalunJ  Pair 

2.7 

Elastic  parameter 

na 

9.66-33.1 

Wing  kinematics 

Formula 

Case  1 

1  -  cos  (<ut) 

Case  2 

sin(6Jt) 

Experimental  data  for  this  case  have  been  received  from  the  University  of  Florida  for  the  10  Hz  case  using 
sin  (cat)  as  the  wing  kinematics.  Grid  sensitivity  analysis  was  conducted  and  the  results  are  shown  in  Figure 
2.6-3,  Based  on  these  results*  the  medium  grid  will  be  employed  for  subsequent  computations. 


Figure  2.6-3.  Time  histories  of  lift  coefficients  (left)  and  vertical  tip  displacement  (right)  for  three  grid  system 
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3.  COMPUTATIONAL  INVESTIGATIONS  OF  FLAPPING  WING 

AERODYNAMICS 


In  this  section,  we  report  the  results  of  the  investigation  of  Reynolds  number  effects  and  shape  effects  on  the 
flow  field  using  two  nominally  two-dimensional  airfoils  (SD7003  and  a  flat  plate)  undergoing  combined 
pitch-  and  plunge  and  pure  plunge  at  Reynolds  numbers  1  *10\  3*I04’  and  6*104.  The  two  different  sets  of 
kinematics  represent  a  weak  dynamic  stall  and  a  stronger  dynamic  stall,  respectively.  Experimental  and 
computational  flowfield  results  are  compared:  phase-averaged  Particle  Image  Velocimetry  (PIV) 
measurements  are  reported,  and  two-dimensional  RANS  equations  coupled  with  Menter’s  Shear  Stress 
Transport  (SST).  In  addition,  lift  coefficient  computed  using  unsteady  linear  airfoil  theory  (Theodorsen 
1935)  is  compared  with  the  computed  lift  coefficient.  The  focus  of  the  investigation  is  to  qualitatively  and 
quantitatively  ascertain  the  role  of  two-dimensional  effects  such  as  leading  edge  vortex  formation,  vortex 
shedding,  and  phase  lag  between  flow  field  and  the  instantaneous  angle  of  attack,  tracing  the  flowfield  and 
lift  coefficient  time  histories.  Issues  such  as  flow  variations  in  the  spanwise  direction,  leading-trailing  edge 
vortex  interaction  with  the  wing  as  well  as  tip  vortices  are  not  addressed  here.  But  as  a  secondary  objective, 
favorable  comparison  between  experiment  and  computation  would  suggest  that  three-dimensional  effects 
would  not  be  of  primary  importance  in  either,  for  the  range  of  motions  presently  under  consideration, 

3. 1  Theodorsen ’s  Unsteady  Linear  Airfoil  Theory 

One  important  issue  in  periodic  oscillatory  airfoil  flows  is  the  lag  between  the  aerodynamic  response  and  the 
airfoil  motion  kinematics.  Quasi-steady  models  for  lift  coefficient  have  enjoyed  some  success  even  in  high- 
frequency  and  geometrical iy-complex  kinematics,  such  as  the  mechanical  models  of  fruit-fly  wings  (Sane 
and  Dickinson  2002).  As  a  natural  extension,  constructing  an  explicit  relation  of  the  lag  of  putatively 
sinusoidal  force  response  to  sinusoidal  motion  kinematics,  as  a  function  of  reduced  frequency,  amplitudes  of 
pitch  and  plunge,  phase  difference  between  pitch  and  plunge,  and  the  Reynolds  number  is  necessary.  This 
could  then  form  a  model  for  the  lift  response  to  more  general  motions  and  in  more  general  configurations. 
Perhaps  the  simplest  generalization  beyond  the  quasi-steady  approximation  was  obtained  by  Theodorsen 
(1935)  model,  for  sinusoidal  pitch-plunge  of  a  thin  airfoil,  by  assuming  a  planar  wake  and  a  trailing-edge 
Kutta  condition,  in  incompressible  inviscid  flow.  The  lift  coefficient  time  history  is  given  by  Eq.  (3,1). 


„  v  ,  .  t ic  (  a  ft  c(2xv  —  ljal 

CM  =  2,(1  -  C(«K  +  y{c+s- 

+  2,C(l.){A+„  +  c(W-2,p)J!-J. 


(3.!) 


The  pitch  and  plunge  motions  are  described  by  the  complex  exponentials,  a(t)  =  a0  +  and 

h(t)  =  h.0e2nf(i.  The  phase  lead  of  pitch  compared  to  plunge  in  terms  of  fractions  of  motion  period  is 
denoted  by  </>.  In  the  most  common  case,  motivated  by  considerations  of  maximum  propulsive  efficiency 
(Anderson  et  al.  1998),  pitch  leads  plunge  by  90°,  which  results  in  (f>  =  0.25.  The  reduced  frequency,  ft,  is 
defined  as  ft  =  nfc/Um  =  nSt/(lhQ),  and  C(ft)  is  the  complex- valued  “Theodorsen  function”  with 
magnitude  <  1.  It  accounts  for  attenuation  of  lift  amplitude  and  time-lag  in  lift  response,  from  its  real  and 
imaginary  parts,  respectively.  The  first  term  is  the  steady-state  lift  and  the  second  term  is  the  "apparent  mass” 
or  noncircuiatory  lift  due  to  acceleration  effects.  The  third  tenn  models  circulatory  effects.  Setting  C(k)  —  1 
(ft  =  0)  recovers  the  quasi-steady  thin  airfoil  solution.  The  noncircuiatory  term  follows  instantaneously  the 
kinematics  of  motion,  but  evolution  of  the  wake  yields  phase  lag  relative  to  the  kinematics  of  airfoil  motion 
in  the  circulatory  tenn,  which  is  predicted  to  peak  for  ft  approximately  equal  to  0.3. 
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The  simplicity  of  Theodorsen’s  model  is  a  powerful  advantage  when  running  targe  parameter  studies,  but  its 
accuracy  for  separated  Hows  with  obviously  nonplanar  wakes  remains  an  issue  of  contention.  In  this  study 
we  compare  the  Theodorsen’s  solution  to  the  RANS  computation  for  lift  coefficient  to  address  the  model’s 
applicability  at  Re  =  0(  104)  for  the  reduced  frequency  of  k  =  0.25. 


3.1.1  Computational  domain  and  kinematics 

The  numerical  solutions  are  computed  in  open  bounded  domain  with  Loci-STREAM  on  an  unstructured  grid 
with  46281,  and  32204  mixed  elements  for  the  SD7003  airfoil,  and  flat  plate,  respectively,  see  Figure  3.1-1* 
The  outer  boundaries  of  the  computational  domain  are  50  (Figure  3.1-l(al)),  and  30  chord  lengths  apart 
(Figure  3.1-1  (bl )),  respectively.  The  thickness  of  the  flat  plate  is  2.3%  chord  length  and  the  leading  and 
trailing  edges  are  rounded  (radius  of  1.15  %  chord  length).  The  boundary  conditions  are  as  follows:  on  the 
airfoil  no-slip  conditions  are  imposed;  the  outer  boundaries  are  incompressible  inlets;  and  the  inlet  turbulence 
intensity  is  0.5%.  The  computations  are  run  assuming  fully-turbulent,  with  no  attempt  to  model  transition  or 
to  prescribe  the  chordwise  location  of  when  to  turn  on  the  production  term  in  the  turbulence  model. 


x/c 

(al)  SD7003  airfoil  in  open  bounded  domain 

mcomp<*»lbtolntet  (MOVING) 


(bl)  Flat  plate  airfoil  in  open  bounded  domain 


(a2)  Mixed  elements  near  the  SD7003  airfoil 


2 


(b2)  Mixed  elements  near  the  flat  plate 


Figure  3.1-1.  Computational  grid  systems:  (a)  SD7003;  (b)  Flat  plate. 


The  motion  kinematics  time  histories  are  described  by 
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h(t)  —  h0c  cos(27t£  /T) 

cr(t)  =  a0  +  Acos{2n(t/T  +  0)} 


where  h  is  the  location  of  the  center  of  rotation  (xp  j c  =  0.25)  of  the  airfoil  measured  normal  to  the  free 
stream,  h0  is  the  normalized  amplitude  of  the  plunge  motion,  T  is  the  motion  physical  period,  c  is  the  airfoil 
chord,  a  is  the  geometrical  angle  of  attack  (AoA)  measured  relative  to  the  incoming  free  stream  with 
velocity,  Um,  o:0  is  the  mean  angle  of  attack,  and  A  is  the  amplitude  of  the  pitching  motion,  see  Figure  3.1-2. 


A  x2 


Pitching 

plunging 


ana 


Figure  3.1-2.  Schematic  of  SD70G3  airfoil  positions  in  downstroke  (RED)  and  upstroke  (BLUE),  and  the 
definition  of  the  free  stream  direction  and  the  effective  angle  of  attack  (effective  AoA)  due  to  plunging  motion. 


The  effective  angle  of  attack,  ae ,  is  a  linear  combination  of  the  pitching  angle  and  the  induced  angle  due  to 
plunging  motion,  and  can  be  written  as, 


ae  =  aQ  +  k  arctan(nSt)  cos{2n(ft  +  0)}  +  arctan{nSt  sin(27r/t)} 


where  St  =  Ifch^/U^  is  the  Strouhal  number,  and  k  —  A/arctan{max(h)/(/m}  is  the  ratio  of  the 
maximum  effective  angles  of  attack  of  the  pitching  motion  to  the  plunge  motion,  where  h  is  the  plunge 
velocity.  Figure  3.1-3.  The  Reynolds  number  is  varied  by  changing  the  flow  speed,  Re  =  V^c/v.  It  is  clear 
from  the  kinematics  that  maintaining  the  same  effective  angle  of  attack  time  history  requires  a  constant 
Strouhal  number  and  constant  A.  Thus,  as  Re  varies,  the  reduced  frequency,  k  =  nfc/U w  —  nStf(2hQ),  and 
the  Strouhal  number  are  kept  constant  by  varying  the  physical  frequency  proportionately. 

The  choice  of  reduced  frequency,  k  =  0.25,  is  motivated  in  part  by  cruise-type  conditions  for  flapping  flight 
of  bird.  Although  the  corresponding  Strouhal  number,  St  =  0.08,  is  below  the  range  for  maximum  propulsion 
efficiency  (Anderson  et  ah  1998),  the  present  flow  conditions  are  on  the  upper-end  of  the  dynamic-stall 
literature,  where  the  main  application  is  helicopter  blade  aerodynamics  (McCroskey  1982)  ,  and  for  which 
the  traditional  analytical  or  phenomenological  models  in  aeronautics  tend  to  focus.  As  is  often  taken  in 
applications  motivated  by  maximizing  propulsive  efficiency  of  pitching  and  plunging  motion,  pitching  leads 
plunging  by  one  quarter  of  motion  period:  phase  0  -  0.25  and  thus  the  airfoil  “feathers"’,  with  the  geometric 
pitching  angle  partially  cancelling  the  plunge-induced  angle  of  attack,  a  rctan  (ft /£/«>).  The  pitching 
amplitude,  A,  is  computed  from  the  value  of  X=  0,6  for  the  combined  pitching  and  plunging  case,  while  for 
the  pure  plunging  ease,  X  =0.  The  total  effective  angle  of  attack  time-trace,  ae,  straddles  the  static  stall  value 
of  approximately  11°  (01  et  aL  2005);  this  is  just  the  sum  of  the  pitching  and  plunging  angles  with 
appropriate  phase  shift. 
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Figure  3.1-3.  Time  history  of  effective  angle  of  attack  (ae)  for  the  pitching  and  plunging  kinematics  (red  line)  and 
the  pure  plunging  kinematics  (blue  line). 


3. 2  Spatial  and  Temporal  Sensitivity  Study 

3.2.1  SD7003 

Spatial  and  temporal  sensitivity  tests  were  performed  for  the  pitch-  and  plunge  case  at  Re=6><104,  k  =  0.25, 
and  0.6,  To  assess  the  grid  sensitivity  time  histories  of  lift  coefficient  on  the  baseline  (46281  cells),  finer 


1  cells), 
itching 


(1 19951  cells)  and  the  finest  (368099  cells)  grids  are  compared  in  Figure  3.2-1  (LEFT)  using  a  time  step  of 
T/dt=  400.  All  three  solutions  coincide,  and  thus  all  subsequent  computations  are  performed  on  the  baseline 
grid.  To  investigate  temporal  sensitivity,  three  time  steps  were  used:  T/dt^  400,  800,  and  1600.  Figure  3.2-1 
(RIGHT)  shows  that  the  computations  using  T/dt  =  400  on  the  grid  with  46281  cells  is  sufficient  to  obtain 
grid  and  time  step-independent  solution. 
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3.2.2  Flat  Plate 

For  the  pitching  and  plunging  flat  plate  the  spatial  sensitivity  test  is  investigated  at  Re  =  6*10\  k  =  0,25,  and 
k  =  0,6,  To  assess  the  grid  sensitivity  time  histories  of  lift  coefficient  on  the  baseline  (9624  cells),  finer 
(32204  cells)  and  the  finest  (65904  cells)  grids  are  compared  in  Figure  3.2-2  using  a  time  step  of  T/dt  =  480. 
All  three  solutions  stay  within  maximum  relative  error  of  2%,  with  the  relative  error  between  the  finer  and 
the  finest  grid  smaller  than  between  the  baseline  and  the  finer  grid.  Based  on  this  observation,  the  finer  grid 
has  been  chosen  for  all  subsequent  computations  for  the  flat  plate. 


Figure  3,2-2,  Time  histories  of  the  lift  coefficients  using  the  baseline  (9624  cells),  finer  (32204  cells),  and  the  finest 
(65904  cells)  grid  using  T/dt  =  480  over  pitching  and  plunging  two-dimensional  flat  plate  at  Re  =  6xlo\k  =  0,25, 
and  k  =  0,6, 


3.3  Flow  around  a  SD7003  Airfoil  at  Re=  6x10* 

3.3.1  Pitching  and  Plunging  Case 

Figure  3.3-1  shows  the  normalized  mean  streamwise  velocity,  Ui/U«,  contours  along  with  planar  streamlines 
from  the  numerical  and  the  experimental  results  from  the  UM  and  AFRL  at  1/T  =  0,00,  0.25,  0,42,  0.50,  and 
0.75,  respectively.  The  numerical  solution  with  the  modified  SST  turbulence  model  overpredicts  the 
separation  leading  to  generation  of  vortical  structures  at  the  bottom  of  the  downstroke,  t/T  =  0.50,  which  is 
not  observed  in  both  PfV  data.  This  is  also  illustrated  in  Figure  3.3-2,  which  shows  U|/U«-component 
velocity  profiles  at  four  different  time  instants  at  constant  x  ,/c=  0.25. 

The  overprediction  of  separation  when  using  the  modified  SST  model  could  be  explained  by  the  use  of  a 
limiter  for  the  production  term  in  the  TKE  equation.  The  build-up  of  turbulence  near  stagnation  flow  region 
is  prevented,  reducing  the  eddy  viscosity  in  the  RANS  model.  Figure  3.3-3  shows  the  local  Reynolds  number 
contours  defined  as  U»c/(v+v,)  from  the  numerical  computations  using  both  SST  turbulence  closures  at  t/T= 
0.25  for  the  pitching  and  plunging  SD7003  airfoil.  The  iimiter  of  the  production  in  the  TKE  equation,  see  Eq. 
(1),  enforced  in  the  modified  SST  model  results  in  substantially  lower  eddy  viscosity,  and  hence  higher  local 
Reynolds  number.  Using  the  original  SST  turbulence  model  the  viscosity  ratio  is  at  maximum  near  the 
leading  edge.  For  the  modified  SST  model,  by  limiting  the  production  of  TKE  the  local  Reynolds  number 
near  the  leading  edge  of  the  airfoil  is  close  to  6*  10\  i.e.  the  amount  of  eddy  viscosity  in  this  region  of  the 
flow  is  small.  Hence  the  flow  tends  to  separate  near  the  leading  edge  which  is  observed  at  t IT-  0.42  and  0.50 
in  Figure  3.3- 1. 

On  the  other  hand,  the  agreement  between  the  two  experimental  measurements  is  excellent,  both  in 
streamwise  velocity  contours  as  well  as  in  streamlines.  During  the  downstroke  motion  the  numerical  solution 
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with  the  modified  SST  model  tends  to  predict  larger  reversed  flow  regions.  The  flow  exhibits  separation 
between  the  center  of  the  downstroke  and  the  bottom  of  the  downstroke  (Figure  3.34),  corresponding  to  the 
maximum  instantaneous  effective  angle  of  attack  of  13.6°.  Note  that  this  value  for  the  effective  angle  of 
attack  is  well  beyond  the  static  stall  angle  of  1 1°. 


Figure  33-1.  u^U*  contours  and  the  instantaneous  streamlines  over  pitching  and  plunging  SD7003  airfoil  at  k  = 
0.25,  k  =0.6,  and  at  Re  =  6*1 04  from  numerical  (Modified  SST,  Original  SST),  and  experimental  (UM,  AFRL) 
results. 
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Figure  33-2*  u^U*  profiles  from  numerical  (Modified  SST,  Original  SST),  and  experimental  (DM,  AFRL) 
results  at  t/T-  0*25,  0*33,  0*42,  and  0.50  at  constant  xi/c=  0*25  at  Re  =  6*104,  k  =  0*25,  X  =0*6  for  the  pitching  and 
plunging  SD7003  airfoil* 


local  Ra:  0  12000  24000  360OD  40000  60000 


(a)  Original  SST  (b)  Modified  SST 


Figure  3*3-3  Local  Reynolds  number  contours  using  (a)  the  original  SST,  and  (b)  the  modified  SST  at  t/T=  0*25 
for  pitching  and  plunging  SD7003  airfoil  at  Re-  6xl04,  k  =  0*25,  X  =0.6. 


Lift  coefficient  time  history  is  shown  in  Figure  3.3-4,  comparing  the  quasi -steady  (2nae%  Theodorsen,  and 
computed  values.  Compared  to  the  steady-state  approximation,  both  Theodorsen *s  result  and  the 
computation  show  smaller  lift  amplitude  as  well  as  some  phase  lag  indicating  non-negligible  influence  from 
the  wake  via  the  circulatory  terms  in  Eq.  (1)  at  k  -  0.25,  Theodorsen’s  solution  and  numerical  solution  agree 
most  closely  at  the  phase  90°.At  phase  180°  the  discrepancy  is  the  largest,  and  hcreboth  numerical  and 
experimental  results  show  an  open  separation  on  the  airfoil  suction  side  (Figure  3*3-1),  Since  the 
Theodorsen5 s  solution  assumes  a  planar  wake  and  Kutta  condition  at  the  trailing  edge,  the  wake  structure  at 
phase  180°  violates  this  condition  causing  the  discrepancy  in  the  lift  coefficient.  Overall,  the  Theodorsen’s 
solution  approximates  the  lift  coefficient  from  the  numerical  compulation  better  when  the  wake  is  'planar’. 
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0.6). 


3.3.2  Pure  plunging  case 


Using  the  original  version  of  SST  turbulence  model  the  computation  showed  a  thinner  but  open  separation, 
however  the  approach  with  the  modified  version  of  SST  model  the  numerical  result  is  able  to  predict  the 
vortical  structure  with  reattachment  at  Xi/c  -  0,8  at  t/T  =  0.25. 


U.AJj  0  0.2  0  4  0.6  0,8  1  1.2 1  4 

Modified  SST 

Original  SST 

UMPIV 

AFRLP1V 

0.00 


Figure  3.3-5.  contours  and  instantaneous  streamlines  over  pure  plunging  S  1)7003  airfoil  at  k=  0.25,  X  =0.0, 
and  at  Re=  6*I0J  from  numerical  (Modified  SST,  Original  SST),  and  experimental  (UM,  AFRL)  results. 
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Figure  3.3-5  shows  the  i^/LL  contour  plots  and  the  instantaneous  streamlines  from  the  numerical 
computation  and  the  experimental  measurements  from  the  UM  and  AFRL  water  tunnels  for  the  pure 
plunging  SD7003  airfoil  at  t/T  -  0.00,  0.25,  0.50,  and  0.75.  The  agreement  between  the  computational  and 
the  experimental  approaches  is  favorable  when  the  flow  is  largely  attached.  When  the  flow  exhibits  massive 
separation,  for  example  at  t/T  =  0.50,  the  experimental  and  computational  results  show  noticeable  differences 
in  phase  as  well  as  the  size  of  flow  separation.  The  details  of  the  vortical  structures  differ  in  ah  results; 
however,  it  is  interesting  to  observe  that  the  original  SST  model  matched  the  PIV  results  from  the  UM 
facility  better,  while  the  modified  SST  model  produced  result  more  consistent  with  that  from  the  AFRL 
facility.  The  consistent/inconsistent  results  appeared  at  t/T  =  0.50  where  a  smaller  vortical  structure  is 
evinced  on  the  suction  side  of  the  airfoil  in  the  UM  facility,  while  in  AFRL  data  such  a  vortical  structure  is 
hardly  present. 


As  already  discussed,  the  flow  tends  to  separate  more  substantially  under  the  modified  SST  model  than  under 
the  original  SST  model  due  to  different  eddy  viscosity  levels  predicted.  The  exact  cause  of  the  difference 
between  the  two  PIV  data  is  not  clear  right  now.  Based  on  the  computational  assessment,  the  effective  inlet 
turbulence  level  of  the  two  tunnels  associated  with  the  wing  motion  may  be  different. 


Figure  33-6  compares  the  lift  coefficient  computed  from  quasi-steady  (2nae%  Theodorsen  and  CFD  for  the 
pure  plunging  case.  Theodorsen's  solution  and  the  numerical  result  coincide  for  t/T=  0.75  to  t/T  =  0.25  while 
between  t/T  =  0.25  and  t/T  -  0.50,  the  numerical  solution  shows  higher  frequency  behavior  and  deviates 
from  the  analytic  prediction  both  in  amplitude  and  phase.  Similar  to  pitching  and  plunging  case,  the  wake 
structures  in  both  PIV  and  CFD  results  are  not  planar  (see  Figure  33-3),  violating  one  of  the  assumptions  for 
the  Theodorsen’s  solution.  The  phase  lag  between  the  effective  angle  of  attack  and  the  response  of  the 
aerodynamic  loading  is  smaller  than  in  the  pitching  and  plunging  case,  despite  the  larger  extent  of  flow 
separation. 


Figure  33-6.  Time  histories  of  lift  coefficient  for  the  pure  plunge  case  (Re=  6><1 04,  k  =  0.25,  X=0.G). 


Unlike  the  pitching  and  plunging  case  where  the  flow  showed  only  thin  open  separation,  the  pure  plunging 
case  generates  large  vortical  structures  at  the  leading  edge  between  motion  phases  of  90°  and  120°. 
Subsequently,  this  structure  -  which  may  be  called  a  leading  edge  vortex  -  broadens,  weakens,  and  convects 
downstream,  eventually  enveloping  the  entire  airfoil  suction  side.  By  180°  phase  of  motion,  reattachment  is 
evinced  at  the  leading  edge,  and  sweeps  downstream  as  the  airfoil  proceeds  on  the  upstroke.  The  LEV  and 
its  subsequent  development  enhance  suction,  and  thus  also  lift.  This  is  seen  in  Figure  33-6  as  a  broad  peak 
in  lift  at  phase  between  90°  and  120°  in  the  numerical  lift  coefficient  result,  followed  by  a  drop  in  lift.  The 
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latter  is  associable  with  weakening  and  downstream  convection  of  the  LEV,  and  loss  of  leading-edge  suction. 
Figure  3.3-7  shows  the  computed  pressure  coefficient  contours,  and  normalized  vorticity  contours  from  both 
the  numerical  and  the  experimental  results  at  the  phase  90°, 


Figure  3*3-7.  Pressure  coefficient  contours  from  the  numerical  computation  (TOP),  normalized  vorticity 
contours  from  the  numerical  computation  (MIDDLE),  and  normalized  vorticity  contours  from  the 
experimental  measurements  (BOTTOM)  at  the  phases  90°  (LEFT),  and  180°  (RIGHT)  for  the  pure  plunge 
case  (Re  =  6*  1 0\  k  -  0.25,  3l=0.0). 


The  LEV  is  notable  in  the  experimental  result,  and  to  a  lesser  extent  in  the  computation.  At  180°  the 
attenuation  in  vorticity  peak  values  is  consistent  with  the  velocity  contour  plots  and  with  the  loss  of  suction 
near  the  leading  edge,  but  there  is  a  notable  discrepancy  between  experiment  and  computation:  the  latter 
shows  a  strong  trailing  edge  vortex,  while  the  former  does  not.  Most  likely,  this  is  the  results  of  poor 
repeatability  of  the  TEV  from  period  to  period,  and  thus  its  dissipation  in  the  phase-averaged  P1V  results. 
Curiously,  the  experimental  and  the  computational  disagreements  seem  to  be  localized  to  the  trailing  edge, 
whence  it  may  be  inferred  that  discrepancy  in  overall  lift  would  be  small  in  the  integrated  sense.  This, 
however,  would  require  substantiation  when  direct  measurement  of  lift  becomes  available  in  the  experiment. 
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